source: trunk/cppad/local/asinh_op.hpp @ 3680

Last change on this file since 3680 was 3680, checked in by bradbell, 5 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 071875a4beba3363e5fa9752426aec4762cd1caa
end hash code: 0bef506513a519e1073c6279d5c4cba9e5c3b180

commit 0bef506513a519e1073c6279d5c4cba9e5c3b180
Author: Brad Bell <bradbell@…>
Date: Thu May 7 12:14:32 2015 -0700

Add the acosh function (as an atomic operation when defined by compiler).

commit b3264fa17b2f65b65800423a0e243c9c3ccfe06a
Author: Brad Bell <bradbell@…>
Date: Wed May 6 20:25:38 2015 -0700

CMakeLists.txt: Change so test only check for compliation.

commit dcbac4d4f20cc383f2bd9edb02036659df40b791
Author: Brad Bell <bradbell@…>
Date: Wed May 6 15:06:28 2015 -0700

asinh.cpp: check higher orders, relax accuracy on test.

commit 5f8881993fedd18cccc3c74831133a8f8a9d17b0
Author: Brad Bell <bradbell@…>
Date: Wed May 6 14:36:18 2015 -0700

Change Acos to acos.
acos.cpp: remove trailing white space.

commit e828fa1f7c4c3848c727f14b1b7a8030071ee705
Author: Brad Bell <bradbell@…>
Date: Wed May 6 12:07:35 2015 -0700

Change Acos to acos.
acos.cpp: remove redundant index commands, remove trailing with space.

commit 3d16e5b9fe1bdafa4ad01d1d466bb72b792650fa
Author: Brad Bell <bradbell@…>
Date: Wed May 6 11:30:49 2015 -0700

op_code.hpp: Minor edits to AcosOp? commnets.

commit 58beaaad149b4ac29fae44589d7f8900bf8f4c40
Author: Brad Bell <bradbell@…>
Date: Wed May 6 10:51:43 2015 -0700

for_jac_sweep.hpp: Add missing AsinhOp? case.

commit 623c134870c522ae5e80bcf0f89d230902594c80
Author: Brad Bell <bradbell@…>
Date: Wed May 6 10:27:39 2015 -0700

Fix comment about AsinhOp? operator.

commit 226b14f6f4810f5abf1ca247aae541963efaf4d6
Author: Brad Bell <bradbell@…>
Date: Wed May 6 10:14:08 2015 -0700

Add derivative of F to make order zero case clearer.
acos_reverse.omh: Fix some sign errors.
asin_reverse.omh: Fix typo.
acos_forward.omh: Simplify by distributing minus sign.

commit 4682f4ee73e33b600b180086576e986f636a24dc
Author: Brad Bell <bradbell@…>
Date: Wed May 6 08:15:50 2015 -0700

acos_forward.omh: fix sign that depends on acos versus acosh.

commit 906ae10adf019ddda7f57dd165aab08fc55289c4
Author: Brad Bell <bradbell@…>
Date: Wed May 6 07:09:47 2015 -0700

  1. Fix inclusion of some temporary files in package (e.g., git_commit.sh).
  2. Simplify and improve using git ls-files and ls bin/check_*.
  3. Remove trailing white space.

commit 5096f4706a547bd76caa3766aa2c62802ef7f0bf
Author: Brad Bell <bradbell@…>
Date: Wed May 6 06:41:20 2015 -0700

Combine base type documentation for erf, asinh
(will add more functions to this list list).

commit b3535db5ad95bee90672abcaa686032d23bce2fc
Author: Brad Bell <bradbell@…>
Date: Tue May 5 18:01:11 2015 -0700

  1. Change Arc Cosine/Sine? to Inverse Cosine/Sine?.
  2. Change arcsin-> asin and arccos->acos.
  3. Remove index commands that are duplicates of words in titles.


acos_reverse.omh: Add acosh case to this page.

File size: 6.8 KB
Line 
1/* $Id$ */
2# ifndef CPPAD_ASINH_OP_INCLUDED
3# define CPPAD_ASINH_OP_INCLUDED
4# if CPPAD_COMPILER_HAS_ASINH
5
6/* --------------------------------------------------------------------------
7CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
8
9CppAD is distributed under multiple licenses. This distribution is under
10the terms of the
11                    Eclipse Public License Version 1.0.
12
13A copy of this license is included in the COPYING file of this distribution.
14Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
15-------------------------------------------------------------------------- */
16
17
18namespace CppAD { // BEGIN_CPPAD_NAMESPACE
19/*!
20\file asinh_op.hpp
21Forward and reverse mode calculations for z = asinh(x).
22*/
23
24
25/*!
26Compute forward mode Taylor coefficient for result of op = AsinhOp.
27
28The C++ source code corresponding to this operation is
29\verbatim
30        z = asinh(x)
31\endverbatim
32The auxillary result is
33\verbatim
34        y = sqrt(1 + x * x)
35\endverbatim
36The value of y, and its derivatives, are computed along with the value
37and derivatives of z.
38
39\copydetails forward_unary2_op
40*/
41template <class Base>
42inline void forward_asinh_op(
43        size_t p           ,
44        size_t q           ,
45        size_t i_z         ,
46        size_t i_x         ,
47        size_t cap_order   ,
48        Base*  taylor      )
49{
50        // check assumptions
51        CPPAD_ASSERT_UNKNOWN( NumArg(AsinhOp) == 1 );
52        CPPAD_ASSERT_UNKNOWN( NumRes(AsinhOp) == 2 );
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        Base uj;
63        if( p == 0 )
64        {       z[0] = asinh( x[0] );
65                uj   = Base(1) + x[0] * x[0];
66                b[0] = sqrt( uj );
67                p++;
68        }
69        for(size_t j = p; j <= q; j++)
70        {       uj = Base(0);
71                for(k = 0; k <= j; k++)
72                        uj += x[k] * x[j-k];
73                b[j] = Base(0);
74                z[j] = Base(0);
75                for(k = 1; k < j; k++)
76                {       b[j] -= Base(k) * b[k] * b[j-k];
77                        z[j] -= Base(k) * z[k] * b[j-k];
78                }
79                b[j] /= Base(j);
80                z[j] /= Base(j);
81                //
82                b[j] += uj / Base(2);
83                z[j] += x[j];
84                //
85                b[j] /= b[0];
86                z[j] /= b[0];
87        }
88}
89/*!
90Multiple directions forward mode Taylor coefficient for op = AsinhOp.
91
92The C++ source code corresponding to this operation is
93\verbatim
94        z = asinh(x)
95\endverbatim
96The auxillary result is
97\verbatim
98        y = sqrt(1 + x * x)
99\endverbatim
100The value of y, and its derivatives, are computed along with the value
101and derivatives of z.
102
103\copydetails forward_unary2_op_dir
104*/
105template <class Base>
106inline void forward_asinh_op_dir(
107        size_t q           ,
108        size_t r           ,
109        size_t i_z         ,
110        size_t i_x         ,
111        size_t cap_order   ,
112        Base*  taylor      )
113{
114        // check assumptions
115        CPPAD_ASSERT_UNKNOWN( NumArg(AcosOp) == 1 );
116        CPPAD_ASSERT_UNKNOWN( NumRes(AcosOp) == 2 );
117        CPPAD_ASSERT_UNKNOWN( 0 < q );
118        CPPAD_ASSERT_UNKNOWN( q < cap_order );
119
120        // Taylor coefficients corresponding to argument and result
121        size_t num_taylor_per_var = (cap_order-1) * r + 1;
122        Base* x = taylor + i_x * num_taylor_per_var;
123        Base* z = taylor + i_z * num_taylor_per_var;
124        Base* b = z - num_taylor_per_var;  // called y in documentation
125
126        size_t k, ell;
127        size_t m = (q-1) * r + 1;
128        for(ell = 0; ell < r; ell ++)
129        {       Base uq = 2.0 * x[m + ell] * x[0];
130                for(k = 1; k < q; k++)
131                        uq += x[(k-1)*r+1+ell] * x[(q-k-1)*r+1+ell];
132                b[m+ell] = Base(0);
133                z[m+ell] = Base(0);
134                for(k = 1; k < q; k++)
135                {       b[m+ell] += Base(k) * b[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell];
136                        z[m+ell] += Base(k) * z[(k-1)*r+1+ell] * b[(q-k-1)*r+1+ell];
137                }
138                b[m+ell] = ( uq / Base(2) - b[m+ell] / Base(q) ) / b[0];
139                z[m+ell] = ( x[m+ell]     - z[m+ell] / Base(q) ) / b[0];
140        }
141}
142
143/*!
144Compute zero order forward mode Taylor coefficient for result of op = AsinhOp.
145
146The C++ source code corresponding to this operation is
147\verbatim
148        z = asinh(x)
149\endverbatim
150The auxillary result is
151\verbatim
152        y = sqrt( 1 + x * x )
153\endverbatim
154The value of y is computed along with the value of z.
155
156\copydetails forward_unary2_op_0
157*/
158template <class Base>
159inline void forward_asinh_op_0(
160        size_t i_z         ,
161        size_t i_x         ,
162        size_t cap_order   ,
163        Base*  taylor      )
164{
165        // check assumptions
166        CPPAD_ASSERT_UNKNOWN( NumArg(AsinhOp) == 1 );
167        CPPAD_ASSERT_UNKNOWN( NumRes(AsinhOp) == 2 );
168        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
169
170        // Taylor coefficients corresponding to argument and result
171        Base* x = taylor + i_x * cap_order;
172        Base* z = taylor + i_z * cap_order;
173        Base* b = z      -       cap_order; // called y in documentation
174
175        z[0] = asinh( x[0] );
176        b[0] = sqrt( Base(1) + x[0] * x[0] );
177}
178/*!
179Compute reverse mode partial derivatives for result of op = AsinhOp.
180
181The C++ source code corresponding to this operation is
182\verbatim
183        z = asinh(x)
184\endverbatim
185The auxillary result is
186\verbatim
187        y = sqrt( 1 + x * x )
188\endverbatim
189The value of y is computed along with the value of z.
190
191\copydetails reverse_unary2_op
192*/
193
194template <class Base>
195inline void reverse_asinh_op(
196        size_t      d            ,
197        size_t      i_z          ,
198        size_t      i_x          ,
199        size_t      cap_order    ,
200        const Base* taylor       ,
201        size_t      nc_partial   ,
202        Base*       partial      )
203{
204        // check assumptions
205        CPPAD_ASSERT_UNKNOWN( NumArg(AsinhOp) == 1 );
206        CPPAD_ASSERT_UNKNOWN( NumRes(AsinhOp) == 2 );
207        CPPAD_ASSERT_UNKNOWN( d < cap_order );
208        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
209
210        // Taylor coefficients and partials corresponding to argument
211        const Base* x  = taylor  + i_x * cap_order;
212        Base* px       = partial + i_x * nc_partial;
213
214        // Taylor coefficients and partials corresponding to first result
215        const Base* z  = taylor  + i_z * cap_order;
216        Base* pz       = partial + i_z * nc_partial;
217
218        // Taylor coefficients and partials corresponding to auxillary result
219        const Base* b  = z  - cap_order; // called y in documentation
220        Base* pb       = pz - nc_partial;
221
222        // If pz is zero, make sure this operation has no effect
223        // (zero times infinity or nan would be non-zero).
224        bool skip(true);
225        for(size_t i_d = 0; i_d <= d; i_d++)
226                skip &= IdenticalZero(pz[i_d]);
227        if( skip )
228                return;
229
230        // number of indices to access
231        size_t j = d;
232        size_t k;
233        while(j)
234        {
235                // scale partials w.r.t b[j] by 1 / b[0]
236                pb[j] /= b[0];
237
238                // scale partials w.r.t z[j] by 1 / b[0]
239                pz[j] /= b[0];
240
241                // update partials w.r.t b^0
242                pb[0] -= pz[j] * z[j] + pb[j] * b[j];
243
244                // update partial w.r.t. x^0
245                px[0] += pb[j] * x[j];
246
247                // update partial w.r.t. x^j
248                px[j] += pz[j] + pb[j] * x[0];
249
250                // further scale partial w.r.t. z[j] by 1 / j
251                pz[j] /= Base(j);
252
253                for(k = 1; k < j; k++)
254                {       // update partials w.r.t b^(j-k)
255                        pb[j-k] -= Base(k) * pz[j] * z[k] + pb[j] * b[k];
256
257                        // update partials w.r.t. x^k
258                        px[k]   += pb[j] * x[j-k];
259
260                        // update partials w.r.t. z^k
261                        pz[k]   -= pz[j] * Base(k) * b[j-k];
262                }
263                --j;
264        }
265
266        // j == 0 case
267        px[0] += ( pz[0] + pb[0] * x[0]) / b[0];
268}
269
270} // END_CPPAD_NAMESPACE
271# endif
272# endif
Note: See TracBrowser for help on using the repository browser.