367
by asaladin
using svn revision number in ptools (1) |
1 |
// $Id$
|
260
by asaladin
copyright for derivify.h |
2 |
/****************************************************************************
|
3 |
* Copyright Joaquim R. R. A. Martins, Peter Sturdza *
|
|
4 |
* *
|
|
5 |
* *
|
|
6 |
* This program is free software: you can redistribute it and/or modify *
|
|
7 |
* it under the terms of the GNU General Public License as published by *
|
|
8 |
* the Free Software Foundation, either version 3 of the License, or *
|
|
9 |
* (at your option) any later version. *
|
|
10 |
* *
|
|
11 |
* This program is distributed in the hope that it will be useful, *
|
|
12 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
|
13 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
|
14 |
* GNU General Public License for more details. *
|
|
15 |
* *
|
|
16 |
* You should have received a copy of the GNU General Public License *
|
|
17 |
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
|
|
18 |
* *
|
|
19 |
***************************************************************************/
|
|
20 |
||
21 |
||
22 |
||
23 |
||
24 |
||
191
by asaladin
informations about the author of derivify.h |
25 |
// NOTE: taken from website http://mdolab.utias.utoronto.ca/resources/complex-step/cpp
|
26 |
//
|
|
260
by asaladin
copyright for derivify.h |
27 |
// please cite: @article{Martins:2003:CSD, Author = {Joaquim R. R. A. Martins and Peter Sturdza and Juan J. Alonso}, Journal = {{ACM} Transactions on Mathematical Software}, Number = {3}, Pages = {245--262}, Title = {The Complex-Step Derivative Approximation}, Volume = {29}, Year = {2003}}
|
28 |
||
29 |
||
190
by asaladin
automatic differentiation header |
30 |
/***
|
31 |
* For automatic differentiation of a C/C++ code:
|
|
32 |
* f'(x) = imag [ f(surreal(x,1)) ]
|
|
33 |
* Define a double precision class that computes and stores values
|
|
34 |
* and derivatives for each variable and overloads all operators.
|
|
35 |
* Tue Jul 18 22:12:28 PDT 2000
|
|
36 |
***/
|
|
37 |
||
38 |
#ifndef HDRderivify
|
|
39 |
#define HDRderivify
|
|
40 |
#include <iostream> |
|
41 |
#include <math.h> |
|
42 |
||
281
by asaladin
modifications to read mbest1u. Force fields parameters now contain dummy atom types |
43 |
// using namespace std;
|
190
by asaladin
automatic differentiation header |
44 |
|
45 |
||
213
by asaladin
multicopy forcefield |
46 |
// #define AUTO_DIFF //uncomment this line to use automatic differenciation
|
194
by asaladin
compile-time switch between double and surreal mode is now working |
47 |
|
192
by asaladin
modifications to use the 'surreal' class to do automatic differentiation. |
48 |
#ifndef HDRcomplexify
|
194
by asaladin
compile-time switch between double and surreal mode is now working |
49 |
inline const double & real(const double& r) { |
233
by asaladin
'astyle' code cleaning applied to all source files |
50 |
/***
|
51 |
* So the real() statement can be used even with
|
|
52 |
* the double version of the code to be derivified,
|
|
53 |
* and remains compatible with complexified code, too.
|
|
54 |
* Most useful inside printf statements.
|
|
55 |
***/
|
|
56 |
return r; |
|
192
by asaladin
modifications to use the 'surreal' class to do automatic differentiation. |
57 |
}
|
58 |
||
194
by asaladin
compile-time switch between double and surreal mode is now working |
59 |
|
60 |
inline double & real(double & r) { return r;} |
|
61 |
||
62 |
||
63 |
||
192
by asaladin
modifications to use the 'surreal' class to do automatic differentiation. |
64 |
inline double imag(const double& r) { |
233
by asaladin
'astyle' code cleaning applied to all source files |
65 |
return 0.; |
192
by asaladin
modifications to use the 'surreal' class to do automatic differentiation. |
66 |
}
|
67 |
#endif // HDRcomplexify |
|
190
by asaladin
automatic differentiation header |
68 |
|
194
by asaladin
compile-time switch between double and surreal mode is now working |
69 |
|
70 |
#ifdef AUTO_DIFF
|
|
71 |
||
190
by asaladin
automatic differentiation header |
72 |
class surreal { |
73 |
#define surr_TEENY (1.e-24) /*machine zero compared to nominal value of val*/ |
|
233
by asaladin
'astyle' code cleaning applied to all source files |
74 |
double val, deriv; |
190
by asaladin
automatic differentiation header |
75 |
|
76 |
public: |
|
233
by asaladin
'astyle' code cleaning applied to all source files |
77 |
surreal(const double& v=0.0, const double& d=0.0) : val(v), deriv(d) {} |
78 |
inline surreal& operator=(const surreal & s) {val = s.val ; deriv = s.deriv; return *this;}; |
|
79 |
operator double() const {return val;} |
|
80 |
operator int() const {return int(val);} |
|
193
by asaladin
small addition to derivify for better python integration |
81 |
|
233
by asaladin
'astyle' code cleaning applied to all source files |
82 |
inline friend const double & real(const surreal& z) ; // born out of |
83 |
inline friend const double & imag(const surreal& z) ; // complex step |
|
84 |
inline friend double & real(surreal& z) ; // born out of |
|
85 |
inline friend double & imag(surreal& z) ; // complex step |
|
86 |
// relational operators
|
|
87 |
inline friend bool operator==(const surreal&,const surreal&); |
|
88 |
inline friend bool operator==(const surreal&,const double&); |
|
89 |
inline friend bool operator==(const double&,const surreal&); |
|
90 |
inline friend bool operator!=(const surreal&,const surreal&); |
|
91 |
inline friend bool operator!=(const surreal&,const double&); |
|
92 |
inline friend bool operator!=(const double&,const surreal&); |
|
93 |
inline friend bool operator>(const surreal&,const surreal&); |
|
94 |
inline friend bool operator>(const surreal&,const double&); |
|
95 |
inline friend bool operator>(const double&,const surreal&); |
|
96 |
inline friend bool operator<(const surreal&,const surreal&); |
|
97 |
inline friend bool operator<(const surreal&,const double&); |
|
98 |
inline friend bool operator<(const double&,const surreal&); |
|
99 |
inline friend bool operator>=(const surreal&,const surreal&); |
|
100 |
inline friend bool operator>=(const surreal&,const double&); |
|
101 |
inline friend bool operator>=(const double&,const surreal&); |
|
102 |
inline friend bool operator<=(const surreal&,const surreal&); |
|
103 |
inline friend bool operator<=(const surreal&,const double&); |
|
104 |
inline friend bool operator<=(const double&,const surreal&); |
|
105 |
// basic arithmetic
|
|
106 |
inline surreal operator+() const; |
|
107 |
inline surreal operator+(const surreal&) const; |
|
108 |
inline surreal operator+(const double&) const; |
|
109 |
inline surreal operator+(const int&) const; |
|
110 |
inline surreal& operator+=(const surreal&); |
|
111 |
inline surreal& operator+=(const double&); |
|
112 |
inline surreal& operator+=(const int&); |
|
113 |
inline friend surreal operator+(const double&, const surreal&); |
|
114 |
inline friend surreal operator+(const int&, const surreal&); |
|
115 |
inline surreal operator-() const; |
|
116 |
inline surreal operator-(const surreal&) const; |
|
117 |
inline surreal operator-(const double&) const; |
|
118 |
inline surreal operator-(const int&) const; |
|
119 |
inline surreal& operator-=(const surreal&); |
|
120 |
inline surreal& operator-=(const double&); |
|
121 |
inline surreal& operator-=(const int&); |
|
122 |
inline friend surreal operator-(const double&, const surreal&); |
|
123 |
inline friend surreal operator-(const int&, const surreal&); |
|
124 |
inline surreal operator*(const surreal&) const; |
|
125 |
inline surreal operator*(const double&) const; |
|
126 |
inline surreal operator*(const int&) const; |
|
127 |
inline surreal& operator*=(const surreal&); |
|
128 |
inline surreal& operator*=(const double&); |
|
129 |
inline surreal& operator*=(const int&); |
|
130 |
inline friend surreal operator*(const double&, const surreal&); |
|
131 |
inline friend surreal operator*(const int&, const surreal&); |
|
132 |
inline surreal operator/(const surreal&) const; |
|
133 |
inline surreal operator/(const double&) const; |
|
134 |
inline surreal operator/(const int&) const; |
|
135 |
inline surreal& operator/=(const surreal&); |
|
136 |
inline surreal& operator/=(const double&); |
|
137 |
inline surreal& operator/=(const int&); |
|
138 |
inline friend surreal operator/(const double&, const surreal&); |
|
139 |
inline friend surreal operator/(const int&, const surreal&); |
|
140 |
// from <math.h>
|
|
141 |
// not implemented are ldexp, frexp, modf, and fmod
|
|
142 |
inline friend surreal fabs(const surreal&); |
|
143 |
inline friend surreal sin(const surreal&); |
|
144 |
inline friend surreal sinh(const surreal&); |
|
145 |
inline friend surreal asin(const surreal&); |
|
146 |
inline friend surreal cos(const surreal&); |
|
147 |
inline friend surreal cosh(const surreal&); |
|
148 |
inline friend surreal acos(const surreal&); |
|
149 |
inline friend surreal tan(const surreal&); |
|
150 |
inline friend surreal tanh(const surreal&); |
|
151 |
inline friend surreal atan(const surreal&); |
|
152 |
inline friend surreal atan2(const surreal&, const surreal&); |
|
153 |
inline friend surreal log(const surreal&); |
|
154 |
inline friend surreal log10(const surreal&); |
|
155 |
inline friend surreal sqrt(const surreal&); |
|
156 |
inline friend surreal exp(const surreal&); |
|
157 |
inline friend surreal pow(const surreal&, const surreal&); |
|
158 |
inline friend surreal pow(const surreal&, const double&); |
|
159 |
inline friend surreal pow(const surreal&, const int&); |
|
160 |
inline friend surreal pow(const double&, const surreal&); |
|
161 |
inline friend surreal pow(const int&, const surreal&); |
|
162 |
inline friend surreal ceil(const surreal&); |
|
163 |
inline friend surreal floor(const surreal&); |
|
164 |
// input/output
|
|
165 |
friend istream& operator>>(istream&, surreal&); |
|
166 |
friend ostream& operator<<(ostream&, const surreal&); |
|
190
by asaladin
automatic differentiation header |
167 |
};
|
168 |
||
169 |
inline bool operator==(const surreal& lhs, const surreal& rhs) |
|
170 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
171 |
return lhs.val == rhs.val; |
190
by asaladin
automatic differentiation header |
172 |
}
|
173 |
||
174 |
inline bool operator==(const surreal& lhs, const double& rhs) |
|
175 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
176 |
return lhs.val == rhs; |
190
by asaladin
automatic differentiation header |
177 |
}
|
178 |
||
179 |
inline bool operator==(const double& lhs, const surreal& rhs) |
|
180 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
181 |
return lhs == rhs.val; |
190
by asaladin
automatic differentiation header |
182 |
}
|
183 |
||
184 |
inline bool operator!=(const surreal& lhs, const surreal& rhs) |
|
185 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
186 |
return lhs.val != rhs.val; |
190
by asaladin
automatic differentiation header |
187 |
}
|
188 |
||
189 |
inline bool operator!=(const surreal& lhs, const double& rhs) |
|
190 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
191 |
return lhs.val != rhs; |
190
by asaladin
automatic differentiation header |
192 |
}
|
193 |
||
194 |
inline bool operator!=(const double& lhs, const surreal& rhs) |
|
195 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
196 |
return lhs != rhs.val; |
190
by asaladin
automatic differentiation header |
197 |
}
|
198 |
||
199 |
inline bool operator>(const surreal& lhs, const surreal& rhs) |
|
200 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
201 |
return lhs.val > rhs.val; |
190
by asaladin
automatic differentiation header |
202 |
}
|
203 |
||
204 |
inline bool operator>(const surreal& lhs, const double& rhs) |
|
205 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
206 |
return lhs.val > rhs; |
190
by asaladin
automatic differentiation header |
207 |
}
|
208 |
||
209 |
inline bool operator>(const double& lhs, const surreal& rhs) |
|
210 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
211 |
return lhs > rhs.val; |
190
by asaladin
automatic differentiation header |
212 |
}
|
213 |
||
214 |
inline bool operator<(const surreal& lhs, const surreal& rhs) |
|
215 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
216 |
return lhs.val < rhs.val; |
190
by asaladin
automatic differentiation header |
217 |
}
|
218 |
||
219 |
inline bool operator<(const surreal& lhs, const double& rhs) |
|
220 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
221 |
return lhs.val < rhs; |
190
by asaladin
automatic differentiation header |
222 |
}
|
223 |
||
224 |
inline bool operator<(const double& lhs, const surreal& rhs) |
|
225 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
226 |
return lhs < rhs.val; |
190
by asaladin
automatic differentiation header |
227 |
}
|
228 |
||
229 |
inline bool operator>=(const surreal& lhs, const surreal& rhs) |
|
230 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
231 |
return lhs.val >= rhs.val; |
190
by asaladin
automatic differentiation header |
232 |
}
|
233 |
||
234 |
inline bool operator>=(const surreal& lhs, const double& rhs) |
|
235 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
236 |
return lhs.val >= rhs; |
190
by asaladin
automatic differentiation header |
237 |
}
|
238 |
||
239 |
inline bool operator>=(const double& lhs, const surreal& rhs) |
|
240 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
241 |
return lhs >= rhs.val; |
190
by asaladin
automatic differentiation header |
242 |
}
|
243 |
||
244 |
inline bool operator<=(const surreal& lhs, const surreal& rhs) |
|
245 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
246 |
return lhs.val <= rhs.val; |
190
by asaladin
automatic differentiation header |
247 |
}
|
248 |
||
249 |
inline bool operator<=(const surreal& lhs, const double& rhs) |
|
250 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
251 |
return lhs.val <= rhs; |
190
by asaladin
automatic differentiation header |
252 |
}
|
253 |
||
254 |
inline bool operator<=(const double& lhs, const surreal& rhs) |
|
255 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
256 |
return lhs <= rhs.val; |
190
by asaladin
automatic differentiation header |
257 |
}
|
258 |
||
259 |
inline surreal surreal::operator+() const |
|
260 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
261 |
return *this; |
190
by asaladin
automatic differentiation header |
262 |
}
|
263 |
||
264 |
inline surreal surreal::operator+(const surreal& z) const |
|
265 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
266 |
return surreal(val+z.val,deriv+z.deriv); |
190
by asaladin
automatic differentiation header |
267 |
}
|
268 |
||
269 |
inline surreal surreal::operator+(const double& r) const |
|
270 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
271 |
return surreal(val+r,deriv); |
190
by asaladin
automatic differentiation header |
272 |
}
|
273 |
||
274 |
inline surreal surreal::operator+(const int& i) const |
|
275 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
276 |
return surreal(val+double(i),deriv); |
190
by asaladin
automatic differentiation header |
277 |
}
|
278 |
||
279 |
inline surreal& surreal::operator+=(const surreal& z) |
|
280 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
281 |
val+=z.val; |
282 |
deriv+=z.deriv; |
|
283 |
return *this; |
|
190
by asaladin
automatic differentiation header |
284 |
}
|
285 |
||
286 |
inline surreal& surreal::operator+=(const double& r) |
|
287 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
288 |
val+=r; |
289 |
return *this; |
|
190
by asaladin
automatic differentiation header |
290 |
}
|
291 |
||
292 |
inline surreal& surreal::operator+=(const int& i) |
|
293 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
294 |
val+=double(i); |
295 |
return *this; |
|
190
by asaladin
automatic differentiation header |
296 |
}
|
297 |
||
298 |
inline surreal operator+(const double& r, const surreal& z) |
|
299 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
300 |
return surreal(r+z.val,z.deriv); |
190
by asaladin
automatic differentiation header |
301 |
}
|
302 |
||
303 |
inline surreal operator+(const int& i, const surreal& z) |
|
304 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
305 |
return surreal(double(i)+z.val,z.deriv); |
190
by asaladin
automatic differentiation header |
306 |
}
|
307 |
||
308 |
inline surreal surreal::operator-() const |
|
309 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
310 |
return surreal(-val,-deriv); |
190
by asaladin
automatic differentiation header |
311 |
}
|
312 |
||
313 |
inline surreal surreal::operator-(const surreal& z) const |
|
314 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
315 |
return surreal(val-z.val,deriv-z.deriv); |
190
by asaladin
automatic differentiation header |
316 |
}
|
317 |
||
318 |
inline surreal surreal::operator-(const double& r) const |
|
319 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
320 |
return surreal(val-r,deriv); |
190
by asaladin
automatic differentiation header |
321 |
}
|
322 |
||
323 |
inline surreal surreal::operator-(const int& i) const |
|
324 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
325 |
return surreal(val-double(i),deriv); |
190
by asaladin
automatic differentiation header |
326 |
}
|
327 |
||
328 |
inline surreal& surreal::operator-=(const surreal& z) |
|
329 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
330 |
val-=z.val; |
331 |
deriv-=z.deriv; |
|
332 |
return *this; |
|
190
by asaladin
automatic differentiation header |
333 |
}
|
334 |
||
335 |
inline surreal& surreal::operator-=(const double& r) |
|
336 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
337 |
val-=r; |
338 |
return *this; |
|
190
by asaladin
automatic differentiation header |
339 |
}
|
340 |
||
341 |
inline surreal& surreal::operator-=(const int& i) |
|
342 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
343 |
val-=double(i); |
344 |
return *this; |
|
190
by asaladin
automatic differentiation header |
345 |
}
|
346 |
||
347 |
inline surreal operator-(const double& r, const surreal& z) |
|
348 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
349 |
return surreal(r-z.val,-z.deriv); |
190
by asaladin
automatic differentiation header |
350 |
}
|
351 |
||
352 |
inline surreal operator-(const int& i, const surreal& z) |
|
353 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
354 |
return surreal(double(i)-z.val,-z.deriv); |
190
by asaladin
automatic differentiation header |
355 |
}
|
356 |
||
357 |
inline surreal surreal::operator*(const surreal& z) const |
|
358 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
359 |
return surreal(val*z.val,val*z.deriv+z.val*deriv); |
190
by asaladin
automatic differentiation header |
360 |
}
|
361 |
||
362 |
inline surreal surreal::operator*(const double& r) const |
|
363 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
364 |
return surreal(val*r,deriv*r); |
190
by asaladin
automatic differentiation header |
365 |
}
|
366 |
||
367 |
inline surreal surreal::operator*(const int& i) const |
|
368 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
369 |
return surreal(val*double(i),deriv*double(i)); |
190
by asaladin
automatic differentiation header |
370 |
}
|
371 |
||
372 |
inline surreal& surreal::operator*=(const surreal& z) |
|
373 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
374 |
deriv=val*z.deriv+z.val*deriv; |
375 |
val*=z.val; |
|
376 |
return *this; |
|
190
by asaladin
automatic differentiation header |
377 |
}
|
378 |
||
379 |
inline surreal& surreal::operator*=(const double& r) |
|
380 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
381 |
val*=r; |
382 |
deriv*=r; |
|
383 |
return *this; |
|
190
by asaladin
automatic differentiation header |
384 |
}
|
385 |
||
386 |
inline surreal& surreal::operator*=(const int& i) |
|
387 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
388 |
val*=double(i); |
389 |
deriv*=double(i); |
|
390 |
return *this; |
|
190
by asaladin
automatic differentiation header |
391 |
}
|
392 |
||
393 |
inline surreal operator*(const double& r, const surreal& z) |
|
394 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
395 |
return surreal(r*z.val,r*z.deriv); |
190
by asaladin
automatic differentiation header |
396 |
}
|
397 |
||
398 |
inline surreal operator*(const int& i, const surreal& z) |
|
399 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
400 |
return surreal(double(i)*z.val,double(i)*z.deriv); |
190
by asaladin
automatic differentiation header |
401 |
}
|
402 |
||
403 |
inline surreal surreal::operator/(const surreal& z) const |
|
404 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
405 |
return surreal(val/z.val,(z.val*deriv-val*z.deriv)/(z.val*z.val)); |
190
by asaladin
automatic differentiation header |
406 |
}
|
407 |
||
408 |
inline surreal surreal::operator/(const double& r) const |
|
409 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
410 |
return surreal(val/r,deriv/r); |
190
by asaladin
automatic differentiation header |
411 |
}
|
412 |
||
413 |
inline surreal surreal::operator/(const int& i) const |
|
414 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
415 |
return surreal(val/double(i),deriv/double(i)); |
190
by asaladin
automatic differentiation header |
416 |
}
|
417 |
||
418 |
inline surreal& surreal::operator/=(const surreal& z) |
|
419 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
420 |
deriv=(z.val*deriv-val*z.deriv)/(z.val*z.val); |
421 |
val/=z.val; |
|
422 |
return *this; |
|
190
by asaladin
automatic differentiation header |
423 |
}
|
424 |
||
425 |
inline surreal& surreal::operator/=(const double& r) |
|
426 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
427 |
val/=r; |
428 |
deriv/=r; |
|
429 |
return *this; |
|
190
by asaladin
automatic differentiation header |
430 |
}
|
431 |
||
432 |
inline surreal& surreal::operator/=(const int& i) |
|
433 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
434 |
val/=double(i); |
435 |
val/=double(i); |
|
436 |
return *this; |
|
190
by asaladin
automatic differentiation header |
437 |
}
|
438 |
||
439 |
inline surreal operator/(const double& r, const surreal& z) |
|
440 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
441 |
return surreal(r/z.val,-r*z.deriv/(z.val*z.val)); |
190
by asaladin
automatic differentiation header |
442 |
}
|
443 |
||
444 |
inline surreal operator/(const int& i, const surreal& z) |
|
445 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
446 |
return surreal(double(i)/z.val,-double(i)*z.deriv/(z.val*z.val)); |
190
by asaladin
automatic differentiation header |
447 |
}
|
448 |
||
449 |
inline surreal fabs(const surreal& z) |
|
450 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
451 |
return (z.val<0.0) ? -z:z; |
190
by asaladin
automatic differentiation header |
452 |
}
|
453 |
||
454 |
inline surreal sin(const surreal& z) |
|
455 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
456 |
return surreal(sin(z.val),z.deriv*cos(z.val)); |
190
by asaladin
automatic differentiation header |
457 |
}
|
458 |
||
459 |
inline surreal sinh(const surreal& z) |
|
460 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
461 |
return surreal(sinh(z.val),z.deriv*cosh(z.val)); |
190
by asaladin
automatic differentiation header |
462 |
}
|
463 |
||
464 |
inline surreal asin(const surreal& z) |
|
465 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
466 |
// derivative trouble if z.val = +/- 1.0
|
467 |
return surreal(asin(z.val),z.deriv/sqrt(1.0-z.val*z.val+surr_TEENY)); |
|
190
by asaladin
automatic differentiation header |
468 |
}
|
469 |
||
470 |
inline surreal cos(const surreal& z) |
|
471 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
472 |
return surreal(cos(z.val),-z.deriv*sin(z.val)); |
190
by asaladin
automatic differentiation header |
473 |
}
|
474 |
||
475 |
inline surreal cosh(const surreal& z) |
|
476 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
477 |
return surreal(cosh(z.val),z.deriv*sinh(z.val)); |
190
by asaladin
automatic differentiation header |
478 |
}
|
479 |
||
480 |
inline surreal acos(const surreal& z) |
|
481 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
482 |
// derivative trouble if z.val = +/- 1.0
|
483 |
return surreal(acos(z.val),-z.deriv/sqrt(1.0-z.val*z.val+surr_TEENY)); |
|
190
by asaladin
automatic differentiation header |
484 |
}
|
485 |
||
486 |
inline surreal tan(const surreal& z) |
|
487 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
488 |
double cosv=cos(z.val); |
489 |
return surreal(tan(z.val),z.deriv/(cosv*cosv)); |
|
190
by asaladin
automatic differentiation header |
490 |
}
|
491 |
||
492 |
inline surreal tanh(const surreal& z) |
|
493 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
494 |
double coshv=cosh(z.val); |
495 |
return surreal(tanh(z.val),z.deriv/(coshv*coshv)); |
|
190
by asaladin
automatic differentiation header |
496 |
}
|
497 |
||
498 |
inline surreal atan(const surreal& z) |
|
499 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
500 |
return surreal(atan(z.val),z.deriv/(1.0+z.val*z.val)); |
190
by asaladin
automatic differentiation header |
501 |
}
|
502 |
||
503 |
inline surreal atan2(const surreal& z1, const surreal& z2) |
|
504 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
505 |
return surreal(atan2(z1.val,z2.val), |
506 |
(z2.val*z1.deriv-z1.val*z2.deriv)/(z1.val*z1.val+z2.val*z2.val)); |
|
190
by asaladin
automatic differentiation header |
507 |
}
|
508 |
||
509 |
inline surreal log(const surreal& z) |
|
510 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
511 |
return surreal(log(z.val),z.deriv/z.val); |
190
by asaladin
automatic differentiation header |
512 |
}
|
513 |
||
514 |
inline surreal log10(const surreal& z) |
|
515 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
516 |
return surreal(log10(z.val),z.deriv/(z.val*log(10.))); |
190
by asaladin
automatic differentiation header |
517 |
}
|
518 |
||
519 |
inline surreal sqrt(const surreal& z) |
|
520 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
521 |
// if z.val = 0, then there is trouble with the derivative.
|
522 |
// this may work if nominal z.val values are scaled properly.
|
|
523 |
double sqrtv=sqrt(z.val); |
|
524 |
return surreal(sqrtv,0.5*z.deriv/(sqrtv+surr_TEENY)); |
|
190
by asaladin
automatic differentiation header |
525 |
}
|
526 |
||
527 |
inline surreal exp(const surreal& z) |
|
528 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
529 |
double expv=exp(z.val); |
530 |
return surreal(expv,z.deriv*expv); |
|
190
by asaladin
automatic differentiation header |
531 |
}
|
532 |
||
533 |
inline surreal pow(const surreal& a, const surreal& b) |
|
534 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
535 |
// many sticky points were derivative is undefined or infinite
|
536 |
// badness if 0 <= b.val < 1 and a.val == 0
|
|
537 |
double powab=pow(a.val,b.val); |
|
538 |
return surreal(powab, |
|
539 |
b.val*pow(a.val,b.val-1.)*a.deriv |
|
540 |
+powab*log(a.val)*b.deriv); |
|
190
by asaladin
automatic differentiation header |
541 |
}
|
542 |
||
543 |
inline surreal pow(const surreal& a, const double& b) |
|
544 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
545 |
return surreal(pow(a.val,b),b*pow(a.val,b-1.)*a.deriv); |
190
by asaladin
automatic differentiation header |
546 |
}
|
547 |
||
548 |
inline surreal pow(const surreal& a, const int& b) |
|
549 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
550 |
return surreal(pow(a.val,double(b)), |
551 |
double(b)*pow(a.val,double(b-1))*a.deriv); |
|
190
by asaladin
automatic differentiation header |
552 |
}
|
553 |
||
554 |
inline surreal pow(const double& a, const surreal& b) |
|
555 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
556 |
double powab=pow(a,b.val); |
557 |
return surreal(powab,powab*log(a)*b.deriv); |
|
190
by asaladin
automatic differentiation header |
558 |
}
|
559 |
||
560 |
inline surreal pow(const int& a, const surreal& b) |
|
561 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
562 |
double powab=pow(double(a),b); |
563 |
return surreal(powab,powab*log(double(a))*b.deriv); |
|
190
by asaladin
automatic differentiation header |
564 |
}
|
565 |
||
566 |
inline surreal ceil(const surreal& z) |
|
567 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
568 |
return surreal(ceil(z.val),0.); |
190
by asaladin
automatic differentiation header |
569 |
}
|
570 |
||
571 |
inline surreal floor(const surreal& z) |
|
572 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
573 |
return surreal(floor(z.val),0.); |
190
by asaladin
automatic differentiation header |
574 |
}
|
575 |
||
576 |
||
577 |
||
578 |
inline istream& operator>>(istream& is, surreal& x) |
|
579 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
580 |
/***
|
581 |
* Straight from Stroustrup's (third edition) complex version.
|
|
582 |
* Much confusion about the ipfx, isfx functions used in
|
|
583 |
* the GNU library and the sentry type that is supposed to
|
|
584 |
* call them indirectly. Stroustrup implies that this
|
|
585 |
* works fine when streams are tied and other fancy stuff.
|
|
586 |
***/
|
|
587 |
||
588 |
double re, im = 0; |
|
589 |
char c = 0; |
|
590 |
||
591 |
is >> c; |
|
592 |
if (c == '(') { |
|
593 |
is >> re >> c; |
|
594 |
if (c == ',') is >> im >> c; |
|
595 |
if (c != ')') is.clear(ios::badbit); |
|
596 |
}
|
|
597 |
else { |
|
598 |
is.putback(c); |
|
599 |
is >> re; |
|
600 |
}
|
|
601 |
||
602 |
if (is) x = surreal (re, im); |
|
603 |
return is; |
|
190
by asaladin
automatic differentiation header |
604 |
}
|
605 |
||
606 |
inline ostream& operator<<(ostream& os, const surreal& z) |
|
607 |
{
|
|
233
by asaladin
'astyle' code cleaning applied to all source files |
608 |
return os << '(' << real (z) << ',' << imag (z) << ')'; |
190
by asaladin
automatic differentiation header |
609 |
}
|
610 |
||
611 |
||
612 |
||
192
by asaladin
modifications to use the 'surreal' class to do automatic differentiation. |
613 |
inline const double & real(const surreal& z) {return z.val;} // born out of |
614 |
inline const double & imag(const surreal& z) {return z.deriv;}// complex step |
|
615 |
inline double & real(surreal& z) {return z.val;} // born out of |
|
616 |
inline double & imag(surreal& z) {return z.deriv;}// complex step |
|
190
by asaladin
automatic differentiation header |
617 |
|
618 |
||
194
by asaladin
compile-time switch between double and surreal mode is now working |
619 |
#endif // #ifdef AUTO_DIFF |
620 |
||
621 |
||
190
by asaladin
automatic differentiation header |
622 |
#undef surr_TEENY
|
623 |
#endif
|
|
624 |