source: trunk/ADOL-C/include/adolc/adouble.h @ 360

Last change on this file since 360 was 360, checked in by kulshres, 7 years ago

implement traced logical operators for condassign and advector

this will help to handle code branching as long as condassign or
advector based generalized conditional assignment is used.

These operators are only used if both operands are badoubles. If the
old behaviour is required then at least one of the two operands must
be a double.

For if-then-else branches the old behaviour should continue and we can do
nothing.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

  • Property svn:keywords set to Author Date Id Revision
File size: 44.3 KB
Line 
1/* ---------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3
4 Revision: $Id: adouble.h 360 2012-10-19 09:22:43Z kulshres $
5 Contents: adouble.h contains the basis for the class of adouble
6           included here are all the possible functions defined on
7           the adouble class.  Notice that, as opposed to ealier versions,
8           both the class adub and the class adouble are derived from a base
9           class (badouble).  See below for further explanation.
10
11 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
12               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel,
13               Kshitij Kulshreshtha
14 
15 This file is part of ADOL-C. This software is provided as open source.
16 Any use, reproduction, or distribution of the software constitutes
17 recipient's acceptance of the terms of the accompanying license file.
18 
19---------------------------------------------------------------------------*/
20
21#if !defined(ADOLC_ADOUBLE_H)
22#define ADOLC_ADOUBLE_H 1
23
24/****************************************************************************/
25/*                                                         THIS FILE IS C++ */
26#ifdef __cplusplus
27
28#include <cstdio>
29#include <cstdlib>
30#include <iostream>
31#include <cmath>
32using std::cout;
33using std::cin;
34using std::cerr;
35using std::ostream;
36using std::istream;
37
38#include <adolc/common.h>
39
40/* NOTICE: There are automatic includes at the end of this file! */
41
42#undef TAPELESS
43#undef SAFE
44#if defined(ADOLC_TAPELESS)
45#  define TAPELESS
46#  undef SAFE
47#endif
48
49#if defined(SAFE_ADOLC_TAPELESS)
50#  define TAPELESS
51#  define SAFE
52#endif
53
54#if !defined(TAPELESS)
55
56/****************************************************************************/
57/*                                             FORWARD DECLARATIONS (TAPES) */
58
59/*--------------------------------------------------------------------------*/
60class adouble;
61class adub;
62class badouble;
63
64/*--------------------------------------------------------------------------*/
65void ADOLC_DLL_EXPORT condassign( double &res, const double &cond,
66                                  const double &arg1, const double &arg2 );
67void ADOLC_DLL_EXPORT condassign( double &res, const double &cond,
68                                  const double &arg );
69
70
71/****************************************************************************/
72/*                                                           CLASS BADOUBLE */
73
74/**
75   The class badouble contains the basic definitions for
76   the arithmetic operations, comparisons, etc.
77   This is a basic class from which the adub and adouble are
78   derived.  Notice that the constructors/destructors for
79   the class badouble are of the trivial variety.  This is the
80   main difference among badoubles, adubs, and adoubles.
81*/
82class ADOLC_DLL_EXPORT badouble {
83protected:
84    locint location;
85    badouble( void ) {};
86    // must be public when using gcc >= 3.4 ( problems with value() )
87    // (see GCC 3.4 Release Series - Changes, New Features, and Fixes)
88    //
89    // badouble( const badouble& a ) {location = a.location;};
90    explicit badouble( locint lo ) {
91        location = lo;
92    };
93
94public:
95    /*--------------------------------------------------------------------------*/
96    badouble( const badouble& a ) {
97        location = a.location;
98    }
99    ;           /* ctor */
100
101    inline locint loc( void ) const;                         /* Helpful stuff */
102
103    /*------------------------------------------------------------------------*/
104    badouble& operator >>= ( double& );                        /* Assignments */
105    badouble& operator <<= ( double );
106    void declareIndependent ();
107    void declareDependent ();
108    badouble& operator = ( double );
109    badouble& operator = ( const badouble& );
110    badouble& operator = ( const adub& );
111    double getValue() const;
112    inline double value() const {
113        return getValue();
114    }
115    void setValue ( const double );
116    /* badouble& operator = ( const adouble& );
117       !!! olvo 991210: was the same as badouble-assignment */
118
119    /*--------------------------------------------------------------------------*/
120    friend ADOLC_DLL_EXPORT std::ostream& operator << ( std::ostream&, const badouble& );  /* IO friends */
121    friend ADOLC_DLL_EXPORT std::istream& operator >> ( std::istream&, const badouble& );
122
123    /*------------------------------------------------------------------------*/
124    badouble& operator += ( double );               /* Operation + Assignment */
125    badouble& operator += ( const badouble& );
126    badouble& operator -= ( double y );
127    badouble& operator -= ( const badouble& );
128    badouble& operator *= ( double );
129    badouble& operator *= ( const badouble& );
130    badouble& operator /= ( double );
131    badouble& operator /= ( const badouble& );
132    /* olvo 991122 n2l: new special op_codes */
133    badouble& operator += ( const adub& );
134    badouble& operator -= ( const adub& );
135
136    /*--------------------------------------------------------------------------*/
137    /* Comparison (friends) */
138    friend ADOLC_DLL_EXPORT adub operator != ( const badouble&, const badouble& );
139    //inline friend int operator != ( const badouble&, const badouble& );
140    inline friend int operator != ( double, const badouble& );
141    friend ADOLC_DLL_EXPORT int operator != ( const badouble&, double );
142    friend ADOLC_DLL_EXPORT adub operator == ( const badouble&, const badouble& );
143    //inline friend int operator == ( const badouble&, const badouble& );
144    inline friend int operator == ( double, const badouble& );
145    friend ADOLC_DLL_EXPORT int operator == ( const badouble&, double );
146    friend ADOLC_DLL_EXPORT adub operator <= ( const badouble&, const badouble& );
147    //inline friend int operator <= ( const badouble&, const badouble& );
148    inline friend int operator <= ( double, const badouble& );
149    friend ADOLC_DLL_EXPORT int operator <= ( const badouble&, double );
150    friend ADOLC_DLL_EXPORT adub operator >= ( const badouble&, const badouble& );
151    //inline friend int operator >= ( const badouble&, const badouble& );
152    inline friend int operator >= ( double, const badouble& );
153    friend ADOLC_DLL_EXPORT int operator >= ( const badouble&, double );
154    friend ADOLC_DLL_EXPORT adub operator >  ( const badouble&, const badouble& );
155    //inline friend int operator >  ( const badouble&, const badouble& );
156    inline friend int operator >  ( double, const badouble& );
157    friend ADOLC_DLL_EXPORT int operator >  ( const badouble&, double );
158    friend ADOLC_DLL_EXPORT adub operator <  ( const badouble&, const badouble& );
159    //inline friend int operator <  ( const badouble&, const badouble& );
160    inline friend int operator <  ( double, const badouble& );
161    friend ADOLC_DLL_EXPORT int operator <  ( const badouble&, double );
162
163
164    /*--------------------------------------------------------------------------*/
165    /* sign operators (friends) */
166    friend ADOLC_DLL_EXPORT adub operator + ( const badouble& x );
167    friend ADOLC_DLL_EXPORT adub operator - ( const badouble& x );
168
169    /*--------------------------------------------------------------------------*/
170    /* binary operators (friends) */
171    friend ADOLC_DLL_EXPORT adub operator + ( const badouble&, const badouble& );
172    friend ADOLC_DLL_EXPORT adub operator + ( double, const badouble& );
173    friend ADOLC_DLL_EXPORT adub operator + ( const badouble&, double );
174    friend ADOLC_DLL_EXPORT adub operator - ( const badouble&, const badouble& );
175    inline friend adub operator - ( const badouble&, double );
176    friend ADOLC_DLL_EXPORT adub operator - ( double, const badouble& );
177    friend ADOLC_DLL_EXPORT adub operator * ( const badouble&, const badouble& );
178    friend ADOLC_DLL_EXPORT adub operator * ( double, const badouble& );
179    inline friend adub operator * ( const badouble&, double );
180    inline friend adub operator / ( const badouble&, double );
181    friend ADOLC_DLL_EXPORT adub operator / ( const badouble&, const badouble& );
182    friend ADOLC_DLL_EXPORT adub operator / ( double, const badouble& );
183
184    /*--------------------------------------------------------------------------*/
185    /* unary operators (friends) */
186    friend ADOLC_DLL_EXPORT adub exp  ( const badouble& );
187    friend ADOLC_DLL_EXPORT adub log  ( const badouble& );
188    friend ADOLC_DLL_EXPORT adub sqrt ( const badouble& );
189    friend ADOLC_DLL_EXPORT adub sin  ( const badouble& );
190    friend ADOLC_DLL_EXPORT adub cos  ( const badouble& );
191    friend ADOLC_DLL_EXPORT adub tan  ( const badouble& );
192    friend ADOLC_DLL_EXPORT adub asin ( const badouble& );
193    friend ADOLC_DLL_EXPORT adub acos ( const badouble& );
194    friend ADOLC_DLL_EXPORT adub atan ( const badouble& );
195
196    /*--------------------------------------------------------------------------*/
197    /* special operators (friends) */
198    friend ADOLC_DLL_EXPORT adouble atan2 ( const badouble&, const badouble& );
199    /* no internal use of condassign: */
200    friend ADOLC_DLL_EXPORT adub    pow   ( const badouble&, double );
201    /* uses condassign internally */
202    friend ADOLC_DLL_EXPORT adouble pow   ( const badouble&, const badouble& );
203    friend ADOLC_DLL_EXPORT adouble pow   ( double, const badouble& );
204    friend ADOLC_DLL_EXPORT adub    log10 ( const badouble& );
205    /* User defined version of logarithm to test extend_quad macro */
206    friend ADOLC_DLL_EXPORT adouble myquad( const badouble& );
207
208    /*--------------------------------------------------------------------------*/
209    /* Additional ANSI C standard Math functions Added by DWJ on 8/6/90 */
210    friend ADOLC_DLL_EXPORT adub sinh  ( const badouble& );
211    friend ADOLC_DLL_EXPORT adub cosh  ( const badouble& );
212    friend ADOLC_DLL_EXPORT adub tanh  ( const badouble& );
213#if defined(ATRIG_ERF)
214    friend ADOLC_DLL_EXPORT adub asinh ( const badouble& );
215    friend ADOLC_DLL_EXPORT adub acosh ( const badouble& );
216    friend ADOLC_DLL_EXPORT adub atanh ( const badouble& );
217#endif
218
219    friend ADOLC_DLL_EXPORT adub fabs  ( const badouble& );
220    friend ADOLC_DLL_EXPORT adub ceil  ( const badouble& );
221    friend ADOLC_DLL_EXPORT adub floor ( const badouble& );
222
223    friend ADOLC_DLL_EXPORT adub fmax ( const badouble&, const badouble& );
224    friend ADOLC_DLL_EXPORT adub fmax ( double, const badouble& );
225    friend ADOLC_DLL_EXPORT adub fmax ( const badouble&, double );
226    friend ADOLC_DLL_EXPORT adub fmin ( const badouble&, const badouble& );
227    friend ADOLC_DLL_EXPORT adub fmin ( double, const badouble& );
228    friend ADOLC_DLL_EXPORT adub fmin ( const badouble&, double );
229
230    friend ADOLC_DLL_EXPORT adub ldexp ( const badouble&, int );
231    friend ADOLC_DLL_EXPORT adub frexp ( const badouble&, int* );
232    friend ADOLC_DLL_EXPORT adub erf   ( const badouble& );
233
234    /*--------------------------------------------------------------------------*/
235    /* Conditionals */
236    friend ADOLC_DLL_EXPORT void condassign( adouble &res, const badouble &cond,
237            const badouble &arg1, const badouble &arg2 );
238    friend ADOLC_DLL_EXPORT void condassign( adouble &res, const badouble &cond,
239            const badouble &arg );
240};
241
242
243
244/****************************************************************************/
245/*                                                               CLASS ADUB */
246
247/*
248   The class Adub
249   ---- Basically used as a temporary result.  The address for an
250        adub is usually generated within an operation.  That address
251        is "freed" when the adub goes out of scope (at destruction time).
252   ---- operates just like a badouble, but it has a destructor defined for it.
253*/
254
255class ADOLC_DLL_EXPORT adub:public badouble {
256    friend ADOLC_DLL_EXPORT class adouble;
257    friend ADOLC_DLL_EXPORT class advector;
258    friend ADOLC_DLL_EXPORT class adubref;
259#if GCC_VERSION >= 4003
260    adub( adub const &) {}
261#endif
262protected:
263    adub( locint lo ):badouble(lo) {};
264    adub( void ):badouble(0) {
265        fprintf(DIAG_OUT,"ADOL-C error: illegal default construction of adub"
266                " variable\n");
267        exit(-2);
268    };
269    explicit adub( double ):badouble(0) {
270        fprintf(DIAG_OUT,"ADOL-C error: illegal  construction of adub variable"
271                " from double\n");
272        exit(-2);
273    };
274
275public:
276
277    /*--------------------------------------------------------------------------*/
278    /* Comparison (friends) */
279    friend ADOLC_DLL_EXPORT adub operator != ( const badouble&, const badouble& );
280    friend ADOLC_DLL_EXPORT adub operator == ( const badouble&, const badouble& );
281    friend ADOLC_DLL_EXPORT adub operator <= ( const badouble&, const badouble& );
282    friend ADOLC_DLL_EXPORT adub operator >= ( const badouble&, const badouble& );
283    friend ADOLC_DLL_EXPORT adub operator < ( const badouble&, const badouble& );
284    friend ADOLC_DLL_EXPORT adub operator > ( const badouble&, const badouble& );
285    /*--------------------------------------------------------------------------*/
286    /* sign operators (friends) */
287    friend ADOLC_DLL_EXPORT adub operator + ( const badouble& x );
288    friend ADOLC_DLL_EXPORT adub operator - ( const badouble& x );
289
290    /*--------------------------------------------------------------------------*/
291    /* binary operators (friends) */
292    friend ADOLC_DLL_EXPORT adub operator + ( const badouble&, const badouble& );
293    friend ADOLC_DLL_EXPORT adub operator + ( double, const badouble& );
294    friend ADOLC_DLL_EXPORT adub operator + ( const badouble&, double );
295    friend ADOLC_DLL_EXPORT adub operator - ( const badouble&, const badouble& );
296    inline friend adub operator - ( const badouble&, double );
297    friend ADOLC_DLL_EXPORT adub operator - ( double, const badouble& );
298    friend ADOLC_DLL_EXPORT adub operator * ( const badouble&, const badouble& );
299    friend ADOLC_DLL_EXPORT adub operator * ( double, const badouble& );
300    inline friend adub operator * ( const badouble&, double );
301    inline friend adub operator / ( const badouble&, double );
302    friend ADOLC_DLL_EXPORT adub operator / ( const badouble&, const badouble& );
303    friend ADOLC_DLL_EXPORT adub operator / ( double, const badouble& );
304
305    /*--------------------------------------------------------------------------*/
306    /* unary operators (friends) */
307    friend ADOLC_DLL_EXPORT adub exp  ( const badouble& );
308    friend ADOLC_DLL_EXPORT adub log  ( const badouble& );
309    friend ADOLC_DLL_EXPORT adub sqrt ( const badouble& );
310    friend ADOLC_DLL_EXPORT adub sin  ( const badouble& );
311    friend ADOLC_DLL_EXPORT adub cos  ( const badouble& );
312    friend ADOLC_DLL_EXPORT adub tan  ( const badouble& );
313    friend ADOLC_DLL_EXPORT adub asin ( const badouble& );
314    friend ADOLC_DLL_EXPORT adub acos ( const badouble& );
315    friend ADOLC_DLL_EXPORT adub atan ( const badouble& );
316
317    /*--------------------------------------------------------------------------*/
318    /* special operators (friends) */
319    /* no internal use of condassign: */
320    friend ADOLC_DLL_EXPORT adub    pow   ( const badouble&, double );
321    friend ADOLC_DLL_EXPORT adub    log10 ( const badouble& );
322
323    /*--------------------------------------------------------------------------*/
324    /* Additional ANSI C standard Math functions Added by DWJ on 8/6/90 */
325    friend ADOLC_DLL_EXPORT adub sinh  ( const badouble& );
326    friend ADOLC_DLL_EXPORT adub cosh  ( const badouble& );
327    friend ADOLC_DLL_EXPORT adub tanh  ( const badouble& );
328#if defined(ATRIG_ERF)
329    friend ADOLC_DLL_EXPORT adub asinh ( const badouble& );
330    friend ADOLC_DLL_EXPORT adub acosh ( const badouble& );
331    friend ADOLC_DLL_EXPORT adub atanh ( const badouble& );
332#endif
333
334    friend ADOLC_DLL_EXPORT adub fabs  ( const badouble& );
335    friend ADOLC_DLL_EXPORT adub ceil  ( const badouble& );
336    friend ADOLC_DLL_EXPORT adub floor ( const badouble& );
337
338    friend ADOLC_DLL_EXPORT adub fmax ( const badouble&, const badouble& );
339    friend ADOLC_DLL_EXPORT adub fmax ( double, const badouble& );
340    friend ADOLC_DLL_EXPORT adub fmax ( const badouble&, double );
341    friend ADOLC_DLL_EXPORT adub fmin ( const badouble&, const badouble& );
342    friend ADOLC_DLL_EXPORT adub fmin ( double, const badouble& );
343    friend ADOLC_DLL_EXPORT adub fmin ( const badouble&, double );
344
345    friend ADOLC_DLL_EXPORT adub ldexp ( const badouble&, int );
346    friend ADOLC_DLL_EXPORT adub frexp ( const badouble&, int* );
347    friend ADOLC_DLL_EXPORT adub erf   ( const badouble& );
348
349    ~adub();
350};
351
352
353/****************************************************************************/
354/*                                                            CLASS ADOUBLE */
355/*
356  The class adouble.
357  ---Derived from badouble.  Contains the standard constructors/destructors.
358  ---At construction, it is given a new address, and at destruction, that
359     address is freed.
360*/
361class ADOLC_DLL_EXPORT adouble:public badouble {
362    friend ADOLC_DLL_EXPORT class advector;
363public:
364    adouble( const adub& );
365    adouble( const adouble& );
366    adouble( void );
367    adouble( double );
368    /* adub prevents postfix operators to occur on the left
369       side of an assignment which would not work  */
370    adub operator++( int );
371    adub operator--( int );
372    badouble& operator++( void );
373    badouble& operator--( void );
374    /*   inline double value(); */
375    ~adouble();
376
377    adouble& operator = ( double );
378    adouble& operator = ( const badouble& );
379    /* adouble& operator = ( const adouble& );
380       !!! olvo 991210 was the same as badouble-assignment */
381    adouble& operator = ( const adub& );
382};
383
384
385/****************************************************************************/
386/*                                                       INLINE DEFINITIONS */
387
388/*--------------------------------------------------------------------------*/
389inline locint badouble::loc( void ) const {
390    return location;
391}
392
393/*--------------------------------------------------------------------------*/
394/* Comparison */
395/*
396inline int operator != ( const badouble& u, const badouble& v ) {
397    return (u-v != 0);
398}
399*/
400
401inline int operator != ( double coval, const badouble& v) {
402    if (coval)
403        return (-coval+v != 0);
404    else
405        return (v != 0);
406}
407
408/*
409inline int operator == ( const badouble& u, const badouble& v ) {
410    return (u-v == 0);
411}
412*/
413
414inline int operator == ( double coval, const badouble& v) {
415    if (coval)
416        return (-coval+v == 0);
417    else
418        return (v == 0);
419}
420
421/*
422inline int operator <= ( const badouble& u, const badouble& v ) {
423    return (u-v <= 0);
424}
425*/
426
427inline int operator <= ( double coval, const badouble& v ) {
428    if (coval)
429        return (-coval+v >= 0);
430    else
431        return (v >= 0);
432}
433
434/*
435inline int operator >= ( const badouble& u, const badouble& v ) {
436    return (u-v >= 0);
437}
438*/
439
440inline int operator >= ( double coval, const badouble& v ) {
441    if (coval)
442        return (-coval+v <= 0);
443    else
444        return (v <= 0);
445}
446
447/*
448inline int operator > ( const badouble& u, const badouble& v ) {
449    return (u-v > 0);
450}
451*/
452
453inline int operator > ( double coval, const badouble& v ) {
454    if (coval)
455        return (-coval+v < 0);
456    else
457        return (v < 0);
458}
459
460/*
461inline int operator < ( const badouble& u, const badouble& v ) {
462    return (u-v < 0);
463}
464*/
465
466inline int operator < ( double coval, const badouble& v ) {
467    if (coval)
468        return (-coval+v > 0);
469    else
470        return (v > 0);
471}
472
473/*--------------------------------------------------------------------------*/
474/* Subtract a floating point from an adouble  */
475inline adub operator - ( const badouble& x , double coval ) {
476    return (-coval) + x;
477}
478
479/*--------------------------------------------------------------------------*/
480/* Multiply an adouble by a floating point */
481inline adub operator * (const badouble& x, double coval) {
482    return coval * x;
483}
484
485/*--------------------------------------------------------------------------*/
486/* Divide an adouble by a floating point */
487inline adub operator / (const badouble& x, double coval) {
488    return (1.0/coval) * x;
489}
490
491/****************************************************************************/
492/* tapeless implementation                                                  */
493/****************************************************************************/
494#else
495
496#include <limits>
497
498namespace adtl {
499
500#if defined(NUMBER_DIRECTIONS)
501extern int ADOLC_numDir;
502#define ADOLC_TAPELESS_UNIQUE_INTERNALS int adtl::ADOLC_numDir = NUMBER_DIRECTIONS;
503#if !defined(DYNAMIC_DIRECTIONS)
504#  define ADVAL                adval[NUMBER_DIRECTIONS]
505#else
506#  define ADVAL                *adval;
507#endif
508#  define ADVAL_TYPE           const double *
509#  define FOR_I_EQ_0_LT_NUMDIR for (int _i=0; _i < ADOLC_numDir; ++_i)
510#  define ADVAL_I              adval[_i]
511#  define ADV_I                adv[_i]
512#  define V_I                  v[_i]
513#else
514#  define ADVAL                adval
515#  define ADVAL_TYPE           double
516#  define FOR_I_EQ_0_LT_NUMDIR
517#  define ADVAL_I              adval
518#  define ADV_I                adv
519#  define V_I                  v
520#endif
521
522#if !defined(_ISOC99_SOURCE) && !defined(__USE_ISOC99)
523inline double fmin( const double &x, const double &y ) {
524    if (x<y) return x;
525    else return y;
526}
527
528inline double fmax( const double &x, const double &y ) {
529    if (x>y) return x;
530    else return y;
531}
532#endif
533
534inline double makeNaN() {
535    return ADOLC_MATH_NSP::numeric_limits<double>::quiet_NaN();
536}
537
538inline double makeInf() {
539    return ADOLC_MATH_NSP::numeric_limits<double>::infinity();
540}
541
542class adouble {
543public:
544    // ctors
545    inline adouble();
546    inline adouble(const double v);
547    inline adouble(const double v, ADVAL_TYPE adv);
548    inline adouble(const adouble& a);
549#if defined(DYNAMIC_DIRECTIONS)
550    inline ~adouble();
551#endif
552
553    /*******************  temporary results  ******************************/
554    // sign
555    inline adouble operator - () const;
556    inline adouble operator + () const;
557
558    // addition
559    inline adouble operator + (const double v) const;
560    inline adouble operator + (const adouble& a) const;
561    inline friend
562    adouble operator + (const double v, const adouble& a);
563
564    // substraction
565    inline adouble operator - (const double v) const;
566    inline adouble operator - (const adouble& a) const;
567    inline friend
568    adouble operator - (const double v, const adouble& a);
569
570    // multiplication
571    inline adouble operator * (const double v) const;
572    inline adouble operator * (const adouble& a) const;
573    inline friend
574    adouble operator * (const double v, const adouble& a);
575
576    // division
577    inline adouble operator / (const double v) const;
578    inline adouble operator / (const adouble& a) const;
579    inline friend
580    adouble operator / (const double v, const adouble& a);
581
582    // inc/dec
583    inline adouble operator ++ ();
584    inline adouble operator ++ (int);
585    inline adouble operator -- ();
586    inline adouble operator -- (int);
587
588    // functions
589    inline friend adouble tan(const adouble &a);
590    inline friend adouble exp(const adouble &a);
591    inline friend adouble log(const adouble &a);
592    inline friend adouble sqrt(const adouble &a);
593    inline friend adouble sin(const adouble &a);
594    inline friend adouble cos(const adouble &a);
595    inline friend adouble asin(const adouble &a);
596    inline friend adouble acos(const adouble &a);
597    inline friend adouble atan(const adouble &a);
598
599    inline friend adouble atan2(const adouble &a, const adouble &b);
600    inline friend adouble pow(const adouble &a, double v);
601    inline friend adouble pow(const adouble &a, const adouble &b);
602    inline friend adouble pow(double v, const adouble &a);
603    inline friend adouble log10(const adouble &a);
604
605    inline friend adouble sinh (const adouble &a);
606    inline friend adouble cosh (const adouble &a);
607    inline friend adouble tanh (const adouble &a);
608#if defined(ATRIG_ERF)
609    inline friend adouble asinh (const adouble &a);
610    inline friend adouble acosh (const adouble &a);
611    inline friend adouble atanh (const adouble &a);
612#endif
613    inline friend adouble fabs (const adouble &a);
614    inline friend adouble ceil (const adouble &a);
615    inline friend adouble floor (const adouble &a);
616    inline friend adouble fmax (const adouble &a, const adouble &b);
617    inline friend adouble fmax (double v, const adouble &a);
618    inline friend adouble fmax (const adouble &a, double v);
619    inline friend adouble fmin (const adouble &a, const adouble &b);
620    inline friend adouble fmin (double v, const adouble &a);
621    inline friend adouble fmin (const adouble &a, double v);
622    inline friend adouble ldexp (const adouble &a, const adouble &b);
623    inline friend adouble ldexp (const adouble &a, const double v);
624    inline friend adouble ldexp (const double v, const adouble &a);
625    inline friend double frexp (const adouble &a, int* v);
626#if defined(ATRIG_ERF)
627    inline friend adouble erf (const adouble &a);
628#endif
629
630
631    /*******************  nontemporary results  ***************************/
632    // assignment
633    inline void operator = (const double v);
634    inline void operator = (const adouble& a);
635
636    // addition
637    inline void operator += (const double v);
638    inline void operator += (const adouble& a);
639
640    // substraction
641    inline void operator -= (const double v);
642    inline void operator -= (const adouble& a);
643
644    // multiplication
645    inline void operator *= (const double v);
646    inline void operator *= (const adouble& a);
647
648    // division
649    inline void operator /= (const double v);
650    inline void operator /= (const adouble& a);
651
652    // not
653    inline int operator ! () const;
654
655    // comparision
656    inline int operator != (const adouble&) const;
657    inline int operator != (const double) const;
658    inline friend int operator != (const double, const adouble&);
659
660    inline int operator == (const adouble&) const;
661    inline int operator == (const double) const;
662    inline friend int operator == (const double, const adouble&);
663
664    inline int operator <= (const adouble&) const;
665    inline int operator <= (const double) const;
666    inline friend int operator <= (const double, const adouble&);
667
668    inline int operator >= (const adouble&) const;
669    inline int operator >= (const double) const;
670    inline friend int operator >= (const double, const adouble&);
671
672    inline int operator >  (const adouble&) const;
673    inline int operator >  (const double) const;
674    inline friend int operator >  (const double, const adouble&);
675
676    inline int operator <  (const adouble&) const;
677    inline int operator <  (const double) const;
678    inline friend int operator <  (const double, const adouble&);
679
680    /*******************  getter / setter  ********************************/
681    inline double getValue() const;
682    inline void setValue(const double v);
683    inline ADVAL_TYPE getADValue() const;
684    inline void setADValue(ADVAL_TYPE v);
685#if defined(NUMBER_DIRECTIONS)
686    inline double getADValue(const unsigned int p) const;
687    inline void setADValue(const unsigned int p, const double v);
688#endif
689
690    /*******************  i/o operations  *********************************/
691    inline friend ostream& operator << ( ostream&, const adouble& );
692    inline friend istream& operator >> ( istream&, adouble& );
693
694
695private:
696    // internal variables
697    double val;
698    double ADVAL;
699};
700
701/*******************************  ctors  ************************************/
702adouble::adouble() {
703#if defined(DYNAMIC_DIRECTIONS)
704    adval = new double[ADOLC_numDir];
705#endif
706}
707
708adouble::adouble(const double v) : val(v) {
709#if defined(DYNAMIC_DIRECTIONS)
710    adval = new double[ADOLC_numDir];
711#endif
712    FOR_I_EQ_0_LT_NUMDIR
713    ADVAL_I = 0.0;
714}
715
716adouble::adouble(const double v, ADVAL_TYPE adv) : val(v) {
717#if defined(DYNAMIC_DIRECTIONS)
718    adval = new double[ADOLC_numDir];
719#endif
720    FOR_I_EQ_0_LT_NUMDIR
721    ADVAL_I=ADV_I;
722}
723
724adouble::adouble(const adouble& a) : val(a.val) {
725#if defined(DYNAMIC_DIRECTIONS)
726    adval = new double[ADOLC_numDir];
727#endif
728    FOR_I_EQ_0_LT_NUMDIR
729    ADVAL_I=a.ADVAL_I;
730}
731
732/*******************************  dtors  ************************************/
733#if defined(DYNAMIC_DIRECTIONS)
734adouble::~adouble() {
735    delete[] adval;
736}
737#endif
738
739/*************************  temporary results  ******************************/
740// sign
741adouble adouble::operator - () const {
742    adouble tmp;
743    tmp.val=-val;
744    FOR_I_EQ_0_LT_NUMDIR
745    tmp.ADVAL_I=-ADVAL_I;
746    return tmp;
747}
748
749adouble adouble::operator + () const {
750    return *this;
751}
752
753// addition
754adouble adouble::operator + (const double v) const {
755    return adouble(val+v, adval);
756}
757
758adouble adouble::operator + (const adouble& a) const {
759    adouble tmp;
760    tmp.val=val+a.val;
761    FOR_I_EQ_0_LT_NUMDIR
762    tmp.ADVAL_I=ADVAL_I+a.ADVAL_I;
763    return tmp;
764}
765
766adouble operator + (const double v, const adouble& a) {
767    return adouble(v+a.val, a.adval);
768}
769
770// subtraction
771adouble adouble::operator - (const double v) const {
772    return adouble(val-v, adval);
773}
774
775adouble adouble::operator - (const adouble& a) const {
776    adouble tmp;
777    tmp.val=val-a.val;
778    FOR_I_EQ_0_LT_NUMDIR
779    tmp.ADVAL_I=ADVAL_I-a.ADVAL_I;
780    return tmp;
781}
782
783adouble operator - (const double v, const adouble& a) {
784    adouble tmp;
785    tmp.val=v-a.val;
786    FOR_I_EQ_0_LT_NUMDIR
787    tmp.ADVAL_I=-a.ADVAL_I;
788    return tmp;
789}
790
791// multiplication
792adouble adouble::operator * (const double v) const {
793    adouble tmp;
794    tmp.val=val*v;
795    FOR_I_EQ_0_LT_NUMDIR
796    tmp.ADVAL_I=ADVAL_I*v;
797    return tmp;
798}
799
800adouble adouble::operator * (const adouble& a) const {
801    adouble tmp;
802    tmp.val=val*a.val;
803    FOR_I_EQ_0_LT_NUMDIR
804    tmp.ADVAL_I=ADVAL_I*a.val+val*a.ADVAL_I;
805    return tmp;
806}
807
808adouble operator * (const double v, const adouble& a) {
809    adouble tmp;
810    tmp.val=v*a.val;
811    FOR_I_EQ_0_LT_NUMDIR
812    tmp.ADVAL_I=v*a.ADVAL_I;
813    return tmp;
814}
815
816// division
817adouble adouble::operator / (const double v) const {
818    adouble tmp;
819    tmp.val=val/v;
820    FOR_I_EQ_0_LT_NUMDIR
821    tmp.ADVAL_I=ADVAL_I/v;
822    return tmp;
823}
824
825adouble adouble::operator / (const adouble& a) const {
826    adouble tmp;
827    tmp.val=val/a.val;
828    FOR_I_EQ_0_LT_NUMDIR
829    tmp.ADVAL_I=(ADVAL_I*a.val-val*a.ADVAL_I)/(a.val*a.val);
830    return tmp;
831}
832
833adouble operator / (const double v, const adouble& a) {
834    adouble tmp;
835    tmp.val=v/a.val;
836    FOR_I_EQ_0_LT_NUMDIR
837    tmp.ADVAL_I=(-v*a.ADVAL_I)/(a.val*a.val);
838    return tmp;
839}
840
841// inc/dec
842adouble adouble::operator ++ () {
843    ++val;
844    return *this;
845}
846
847adouble adouble::operator ++ (int) {
848    adouble tmp;
849    tmp.val=val++;
850    FOR_I_EQ_0_LT_NUMDIR
851    tmp.ADVAL_I=ADVAL_I;
852    return tmp;
853}
854
855adouble adouble::operator -- () {
856    --val;
857    return *this;
858}
859
860adouble adouble::operator -- (int) {
861    adouble tmp;
862    tmp.val=val--;
863    FOR_I_EQ_0_LT_NUMDIR
864    tmp.ADVAL_I=ADVAL_I;
865    return tmp;
866}
867
868// functions
869adouble tan(const adouble& a) {
870    adouble tmp;
871    double tmp2;
872    tmp.val=ADOLC_MATH_NSP::tan(a.val);
873    tmp2=ADOLC_MATH_NSP::cos(a.val);
874    tmp2*=tmp2;
875    FOR_I_EQ_0_LT_NUMDIR
876    tmp.ADVAL_I=a.ADVAL_I/tmp2;
877    return tmp;
878}
879
880adouble exp(const adouble &a) {
881    adouble tmp;
882    tmp.val=ADOLC_MATH_NSP::exp(a.val);
883    FOR_I_EQ_0_LT_NUMDIR
884    tmp.ADVAL_I=tmp.val*a.ADVAL_I;
885    return tmp;
886}
887
888adouble log(const adouble &a) {
889    adouble tmp;
890    tmp.val=ADOLC_MATH_NSP::log(a.val);
891    FOR_I_EQ_0_LT_NUMDIR
892        if (a.val>0) tmp.ADVAL_I=a.ADVAL_I/a.val;
893        else if (a.val==0 && a.ADVAL_I != 0.0) {
894            int sign = (a.ADVAL_I < 0)  ? -1 : 1;
895            tmp.ADVAL_I=sign*makeInf();
896        } else tmp.ADVAL_I=makeNaN();
897    return tmp;
898}
899
900adouble sqrt(const adouble &a) {
901    adouble tmp;
902    tmp.val=ADOLC_MATH_NSP::sqrt(a.val);
903    FOR_I_EQ_0_LT_NUMDIR
904        if (a.val>0) tmp.ADVAL_I=a.ADVAL_I/(tmp.val*2);
905        else if (a.val==0.0 && a.ADVAL_I != 0.0) {
906            int sign = (a.ADVAL_I < 0) ? -1 : 1;
907            tmp.ADVAL_I=sign * makeInf();
908        } else tmp.ADVAL_I=makeNaN();
909    return tmp;
910}
911
912adouble sin(const adouble &a) {
913    adouble tmp;
914    double tmp2;
915    tmp.val=ADOLC_MATH_NSP::sin(a.val);
916    tmp2=ADOLC_MATH_NSP::cos(a.val);
917    FOR_I_EQ_0_LT_NUMDIR
918    tmp.ADVAL_I=tmp2*a.ADVAL_I;
919    return tmp;
920}
921
922adouble cos(const adouble &a) {
923    adouble tmp;
924    double tmp2;
925    tmp.val=ADOLC_MATH_NSP::cos(a.val);
926    tmp2=-ADOLC_MATH_NSP::sin(a.val);
927    FOR_I_EQ_0_LT_NUMDIR
928    tmp.ADVAL_I=tmp2*a.ADVAL_I;
929    return tmp;
930}
931
932adouble asin(const adouble &a) {
933    adouble tmp;
934    tmp.val=ADOLC_MATH_NSP::asin(a.val);
935    double tmp2=ADOLC_MATH_NSP::sqrt(1-a.val*a.val);
936    FOR_I_EQ_0_LT_NUMDIR
937    tmp.ADVAL_I=a.ADVAL_I/tmp2;
938    return tmp;
939}
940
941adouble acos(const adouble &a) {
942    adouble tmp;
943    tmp.val=ADOLC_MATH_NSP::acos(a.val);
944    double tmp2=-ADOLC_MATH_NSP::sqrt(1-a.val*a.val);
945    FOR_I_EQ_0_LT_NUMDIR
946    tmp.ADVAL_I=a.ADVAL_I/tmp2;
947    return tmp;
948}
949
950adouble atan(const adouble &a) {
951    adouble tmp;
952    tmp.val=ADOLC_MATH_NSP::atan(a.val);
953    double tmp2=1+a.val*a.val;
954    tmp2=1/tmp2;
955    if (tmp2!=0)
956        FOR_I_EQ_0_LT_NUMDIR
957        tmp.ADVAL_I=a.ADVAL_I*tmp2;
958    else
959        FOR_I_EQ_0_LT_NUMDIR
960        tmp.ADVAL_I=0.0;
961    return tmp;
962}
963
964adouble atan2(const adouble &a, const adouble &b) {
965    adouble tmp;
966    tmp.val=ADOLC_MATH_NSP::atan2(a.val, b.val);
967    double tmp2=a.val*a.val;
968    double tmp3=b.val*b.val;
969    double tmp4=tmp3/(tmp2+tmp3);
970    if (tmp4!=0)
971        FOR_I_EQ_0_LT_NUMDIR
972        tmp.ADVAL_I=(a.ADVAL_I*b.val-a.val*b.ADVAL_I)/tmp3*tmp4;
973    else
974        FOR_I_EQ_0_LT_NUMDIR
975        tmp.ADVAL_I=0.0;
976    return tmp;
977}
978
979adouble pow(const adouble &a, double v) {
980    adouble tmp;
981    tmp.val=ADOLC_MATH_NSP::pow(a.val, v);
982    double tmp2=v*ADOLC_MATH_NSP::pow(a.val, v-1);
983    FOR_I_EQ_0_LT_NUMDIR
984    tmp.ADVAL_I=tmp2*a.ADVAL_I;
985    return tmp;
986}
987
988adouble pow(const adouble &a, const adouble &b) {
989    adouble tmp;
990    tmp.val=ADOLC_MATH_NSP::pow(a.val, b.val);
991    double tmp2=b.val*ADOLC_MATH_NSP::pow(a.val, b.val-1);
992    double tmp3=ADOLC_MATH_NSP::log(a.val)*tmp.val;
993    FOR_I_EQ_0_LT_NUMDIR
994    tmp.ADVAL_I=tmp2*a.ADVAL_I+tmp3*b.ADVAL_I;
995    return tmp;
996}
997
998adouble pow(double v, const adouble &a) {
999    adouble tmp;
1000    tmp.val=ADOLC_MATH_NSP::pow(v, a.val);
1001    double tmp2=tmp.val*ADOLC_MATH_NSP::log(v);
1002    FOR_I_EQ_0_LT_NUMDIR
1003    tmp.ADVAL_I=tmp2*a.ADVAL_I;
1004    return tmp;
1005}
1006
1007adouble log10(const adouble &a) {
1008    adouble tmp;
1009    tmp.val=ADOLC_MATH_NSP::log10(a.val);
1010    double tmp2=ADOLC_MATH_NSP::log((double)10)*a.val;
1011    FOR_I_EQ_0_LT_NUMDIR
1012    tmp.ADVAL_I=a.ADVAL_I/tmp2;
1013    return tmp;
1014}
1015
1016adouble sinh (const adouble &a) {
1017    adouble tmp;
1018    tmp.val=ADOLC_MATH_NSP::sinh(a.val);
1019    double tmp2=ADOLC_MATH_NSP::cosh(a.val);
1020    FOR_I_EQ_0_LT_NUMDIR
1021    tmp.ADVAL_I=a.ADVAL_I*tmp2;
1022    return tmp;
1023}
1024
1025adouble cosh (const adouble &a) {
1026    adouble tmp;
1027    tmp.val=ADOLC_MATH_NSP::cosh(a.val);
1028    double tmp2=ADOLC_MATH_NSP::sinh(a.val);
1029    FOR_I_EQ_0_LT_NUMDIR
1030    tmp.ADVAL_I=a.ADVAL_I*tmp2;
1031    return tmp;
1032}
1033
1034adouble tanh (const adouble &a) {
1035    adouble tmp;
1036    tmp.val=ADOLC_MATH_NSP::tanh(a.val);
1037    double tmp2=ADOLC_MATH_NSP::cosh(a.val);
1038    tmp2*=tmp2;
1039    FOR_I_EQ_0_LT_NUMDIR
1040    tmp.ADVAL_I=a.ADVAL_I/tmp2;
1041    return tmp;
1042}
1043
1044#if defined(ATRIG_ERF)
1045adouble asinh (const adouble &a) {
1046    adouble tmp;
1047    tmp.val=ADOLC_MATH_NSP_ERF::asinh(a.val);
1048    double tmp2=ADOLC_MATH_NSP::sqrt(a.val*a.val+1);
1049    FOR_I_EQ_0_LT_NUMDIR
1050    tmp.ADVAL_I=a.ADVAL_I/tmp2;
1051    return tmp;
1052}
1053
1054adouble acosh (const adouble &a) {
1055    adouble tmp;
1056    tmp.val=ADOLC_MATH_NSP_ERF::acosh(a.val);
1057    double tmp2=ADOLC_MATH_NSP::sqrt(a.val*a.val-1);
1058    FOR_I_EQ_0_LT_NUMDIR
1059    tmp.ADVAL_I=a.ADVAL_I/tmp2;
1060    return tmp;
1061}
1062
1063adouble atanh (const adouble &a) {
1064    adouble tmp;
1065    tmp.val=ADOLC_MATH_NSP_ERF::atanh(a.val);
1066    double tmp2=1-a.val*a.val;
1067    FOR_I_EQ_0_LT_NUMDIR
1068    tmp.ADVAL_I=a.ADVAL_I/tmp2;
1069    return tmp;
1070}
1071#endif
1072
1073adouble fabs (const adouble &a) {
1074    adouble tmp;
1075    tmp.val=ADOLC_MATH_NSP::fabs(a.val);
1076    int as=0;
1077    if (a.val>0) as=1;
1078    if (a.val<0) as=-1;
1079    if (as!=0)
1080        FOR_I_EQ_0_LT_NUMDIR
1081        tmp.ADVAL_I=a.ADVAL_I*as;
1082    else
1083        FOR_I_EQ_0_LT_NUMDIR {
1084            as=0;
1085            if (a.ADVAL_I>0) as=1;
1086            if (a.ADVAL_I<0) as=-1;
1087                tmp.ADVAL_I=a.ADVAL_I*as;
1088            }
1089            return tmp;
1090}
1091
1092adouble ceil (const adouble &a) {
1093    adouble tmp;
1094    tmp.val=ADOLC_MATH_NSP::ceil(a.val);
1095    FOR_I_EQ_0_LT_NUMDIR
1096    tmp.ADVAL_I=0.0;
1097    return tmp;
1098}
1099
1100adouble floor (const adouble &a) {
1101    adouble tmp;
1102    tmp.val=ADOLC_MATH_NSP::floor(a.val);
1103    FOR_I_EQ_0_LT_NUMDIR
1104    tmp.ADVAL_I=0.0;
1105    return tmp;
1106}
1107
1108adouble fmax (const adouble &a, const adouble &b) {
1109    adouble tmp;
1110    double tmp2=a.val-b.val;
1111    if (tmp2<0) {
1112        tmp.val=b.val;
1113        FOR_I_EQ_0_LT_NUMDIR
1114        tmp.ADVAL_I=b.ADVAL_I;
1115    } else {
1116        tmp.val=a.val;
1117        if (tmp2>0) {
1118            FOR_I_EQ_0_LT_NUMDIR
1119            tmp.ADVAL_I=a.ADVAL_I;
1120        } else {
1121            FOR_I_EQ_0_LT_NUMDIR
1122            {
1123                if (a.ADVAL_I<b.ADVAL_I) tmp.ADVAL_I=b.ADVAL_I;
1124                else tmp.ADVAL_I=a.ADVAL_I;
1125                }
1126            }
1127}
1128return tmp;
1129}
1130
1131adouble fmax (double v, const adouble &a) {
1132    adouble tmp;
1133    double tmp2=v-a.val;
1134    if (tmp2<0) {
1135        tmp.val=a.val;
1136        FOR_I_EQ_0_LT_NUMDIR
1137        tmp.ADVAL_I=a.ADVAL_I;
1138    } else {
1139        tmp.val=v;
1140        if (tmp2>0) {
1141            FOR_I_EQ_0_LT_NUMDIR
1142            tmp.ADVAL_I=0.0;
1143        } else {
1144            FOR_I_EQ_0_LT_NUMDIR
1145            {
1146                if (a.ADVAL_I>0) tmp.ADVAL_I=a.ADVAL_I;
1147                else tmp.ADVAL_I=0.0;
1148                }
1149            }
1150}
1151return tmp;
1152}
1153
1154adouble fmax (const adouble &a, double v) {
1155    adouble tmp;
1156    double tmp2=a.val-v;
1157    if (tmp2<0) {
1158        tmp.val=v;
1159        FOR_I_EQ_0_LT_NUMDIR
1160        tmp.ADVAL_I=0.0;
1161    } else {
1162        tmp.val=a.val;
1163        if (tmp2>0) {
1164            FOR_I_EQ_0_LT_NUMDIR
1165            tmp.ADVAL_I=a.ADVAL_I;
1166        } else {
1167            FOR_I_EQ_0_LT_NUMDIR
1168            {
1169                if (a.ADVAL_I>0) tmp.ADVAL_I=a.ADVAL_I;
1170                else tmp.ADVAL_I=0.0;
1171                }
1172            }
1173}
1174return tmp;
1175}
1176
1177adouble fmin (const adouble &a, const adouble &b) {
1178    adouble tmp;
1179    double tmp2=a.val-b.val;
1180    if (tmp2<0) {
1181        tmp.val=a.val;
1182        FOR_I_EQ_0_LT_NUMDIR
1183        tmp.ADVAL_I=a.ADVAL_I;
1184    } else {
1185        tmp.val=b.val;
1186        if (tmp2>0) {
1187            FOR_I_EQ_0_LT_NUMDIR
1188            tmp.ADVAL_I=b.ADVAL_I;
1189        } else {
1190            FOR_I_EQ_0_LT_NUMDIR
1191            {
1192                if (a.ADVAL_I<b.ADVAL_I) tmp.ADVAL_I=a.ADVAL_I;
1193                else tmp.ADVAL_I=b.ADVAL_I;
1194                }
1195            }
1196}
1197return tmp;
1198}
1199
1200adouble fmin (double v, const adouble &a) {
1201    adouble tmp;
1202    double tmp2=v-a.val;
1203    if (tmp2<0) {
1204        tmp.val=v;
1205        FOR_I_EQ_0_LT_NUMDIR
1206        tmp.ADVAL_I=0.0;
1207    } else {
1208        tmp.val=a.val;
1209        if (tmp2>0) {
1210            FOR_I_EQ_0_LT_NUMDIR
1211            tmp.ADVAL_I=a.ADVAL_I;
1212        } else {
1213            FOR_I_EQ_0_LT_NUMDIR
1214            {
1215                if (a.ADVAL_I<0) tmp.ADVAL_I=a.ADVAL_I;
1216                else tmp.ADVAL_I=0.0;
1217                }
1218            }
1219}
1220return tmp;
1221}
1222
1223adouble fmin (const adouble &a, double v) {
1224    adouble tmp;
1225    double tmp2=a.val-v;
1226    if (tmp2<0) {
1227        tmp.val=a.val;
1228        FOR_I_EQ_0_LT_NUMDIR
1229        tmp.ADVAL_I=a.ADVAL_I;
1230    } else {
1231        tmp.val=v;
1232        if (tmp2>0) {
1233            FOR_I_EQ_0_LT_NUMDIR
1234            tmp.ADVAL_I=0.0;
1235        } else {
1236            FOR_I_EQ_0_LT_NUMDIR
1237            {
1238                if (a.ADVAL_I<0) tmp.ADVAL_I=a.ADVAL_I;
1239                else tmp.ADVAL_I=0.0;
1240                }
1241            }
1242}
1243return tmp;
1244}
1245
1246adouble ldexp (const adouble &a, const adouble &b) {
1247    return a*pow(2.,b);
1248}
1249
1250adouble ldexp (const adouble &a, const double v) {
1251    return a*ADOLC_MATH_NSP::pow(2.,v);
1252}
1253
1254adouble ldexp (const double v, const adouble &a) {
1255    return v*pow(2.,a);
1256}
1257
1258double frexp (const adouble &a, int* v) {
1259    return ADOLC_MATH_NSP::frexp(a.val, v);
1260}
1261
1262#if defined(ATRIG_ERF)
1263adouble erf (const adouble &a) {
1264    adouble tmp;
1265    tmp.val=ADOLC_MATH_NSP_ERF::erf(a.val);
1266    double tmp2 = 2.0 /
1267        ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) *
1268        ADOLC_MATH_NSP_ERF::exp(-a.val*a.val);
1269    FOR_I_EQ_0_LT_NUMDIR
1270    tmp.ADVAL_I=tmp2*a.ADVAL_I;
1271    return tmp;
1272}
1273#endif
1274
1275
1276/*******************  nontemporary results  *********************************/
1277void adouble::operator = (const double v) {
1278    val=v;
1279    FOR_I_EQ_0_LT_NUMDIR
1280    ADVAL_I=0.0;
1281}
1282
1283void adouble::operator = (const adouble& a) {
1284    val=a.val;
1285    FOR_I_EQ_0_LT_NUMDIR
1286    ADVAL_I=a.ADVAL_I;
1287}
1288
1289void adouble::operator += (const double v) {
1290    val+=v;
1291}
1292
1293void adouble::operator += (const adouble& a) {
1294    val=val+a.val;
1295    FOR_I_EQ_0_LT_NUMDIR
1296    ADVAL_I+=a.ADVAL_I;
1297}
1298
1299void adouble::operator -= (const double v) {
1300    val-=v;
1301}
1302
1303void adouble::operator -= (const adouble& a) {
1304    val=val-a.val;
1305    FOR_I_EQ_0_LT_NUMDIR
1306    ADVAL_I-=a.ADVAL_I;
1307}
1308
1309void adouble::operator *= (const double v) {
1310    val=val*v;
1311    FOR_I_EQ_0_LT_NUMDIR
1312    ADVAL_I*=v;
1313}
1314
1315void adouble::operator *= (const adouble& a) {
1316    FOR_I_EQ_0_LT_NUMDIR
1317    ADVAL_I=ADVAL_I*a.val+val*a.ADVAL_I;
1318    val*=a.val;
1319}
1320
1321void adouble::operator /= (const double v) {
1322    val/=v;
1323    FOR_I_EQ_0_LT_NUMDIR
1324    ADVAL_I/=v;
1325}
1326
1327void adouble::operator /= (const adouble& a) {
1328    FOR_I_EQ_0_LT_NUMDIR
1329    ADVAL_I=(ADVAL_I*a.val-val*a.ADVAL_I)/(a.val*a.val);
1330    val=val/a.val;
1331}
1332
1333// not
1334int adouble::operator ! () const {
1335    return val==0.0;
1336}
1337
1338// comparision
1339int adouble::operator != (const adouble &a) const {
1340    return val!=a.val;
1341}
1342
1343int adouble::operator != (const double v) const {
1344    return val!=v;
1345}
1346
1347int operator != (const double v, const adouble &a) {
1348    return v!=a.val;
1349}
1350
1351int adouble::operator == (const adouble &a) const {
1352    return val==a.val;
1353}
1354
1355int adouble::operator == (const double v) const {
1356    return val==v;
1357}
1358
1359int operator == (const double v, const adouble &a) {
1360    return v==a.val;
1361}
1362
1363int adouble::operator <= (const adouble &a) const {
1364    return val<=a.val;
1365}
1366
1367int adouble::operator <= (const double v) const {
1368    return val<=v;
1369}
1370
1371int operator <= (const double v, const adouble &a) {
1372    return v<=a.val;
1373}
1374
1375int adouble::operator >= (const adouble &a) const {
1376    return val>=a.val;
1377}
1378
1379int adouble::operator >= (const double v) const {
1380    return val>=v;
1381}
1382
1383int operator >= (const double v, const adouble &a) {
1384    return v>=a.val;
1385}
1386
1387int adouble::operator >  (const adouble &a) const {
1388    return val>a.val;
1389}
1390
1391int adouble::operator >  (const double v) const {
1392    return val>v;
1393}
1394
1395int operator >  (const double v, const adouble &a) {
1396    return v>a.val;
1397}
1398
1399int adouble::operator <  (const adouble &a) const {
1400    return val<a.val;
1401}
1402
1403int adouble::operator <  (const double v) const {
1404    return val<v;
1405}
1406
1407int operator <  (const double v, const adouble &a) {
1408    return v<a.val;
1409}
1410
1411/*******************  getter / setter  **************************************/
1412double adouble::getValue() const {
1413    return val;
1414}
1415
1416void adouble::setValue(const double v) {
1417    val=v;
1418}
1419
1420ADVAL_TYPE adouble::getADValue() const {
1421    return adval;
1422}
1423
1424void adouble::setADValue(ADVAL_TYPE v) {
1425    FOR_I_EQ_0_LT_NUMDIR
1426    ADVAL_I=V_I;
1427}
1428
1429#  if defined(NUMBER_DIRECTIONS)
1430double adouble::getADValue(const unsigned int p) const {
1431    if (p>=NUMBER_DIRECTIONS) {
1432        fprintf(DIAG_OUT, "Derivative array accessed out of bounds"\
1433                " while \"getADValue(...)\"!!!\n");
1434        exit(-1);
1435    }
1436    return adval[p];
1437}
1438
1439void adouble::setADValue(const unsigned int p, const double v) {
1440    if (p>=NUMBER_DIRECTIONS) {
1441        fprintf(DIAG_OUT, "Derivative array accessed out of bounds"\
1442                " while \"setADValue(...)\"!!!\n");
1443        exit(-1);
1444    }
1445    adval[p]=v;
1446}
1447#  endif
1448
1449#if defined(NUMBER_DIRECTIONS)
1450static void setNumDir(const unsigned int p) {
1451#if !defined(DYNAMIC_DIRECTIONS)
1452    if (p>NUMBER_DIRECTIONS) ADOLC_numDir=NUMBER_DIRECTIONS;
1453    else ADOLC_numDir=p;
1454#else
1455    ADOLC_numDir = p;
1456#endif
1457}
1458#endif
1459
1460/*******************  i/o operations  ***************************************/
1461ostream& operator << ( ostream& out, const adouble& a) {
1462    out << "Value: " << a.val;
1463#if !defined(NUMBER_DIRECTIONS)
1464    out << " ADValue: ";
1465#else
1466    out << " ADValues (" << ADOLC_numDir << "): ";
1467#endif
1468    FOR_I_EQ_0_LT_NUMDIR
1469    out << a.ADVAL_I << " ";
1470    out << "(a)";
1471    return out;
1472}
1473
1474istream& operator >> ( istream& in, adouble& a) {
1475    char c;
1476    do {
1477        in >> c;
1478    } while (c!=':' && !in.eof());
1479    in >> a.val;
1480#if !defined(NUMBER_DIRECTIONS)
1481    do in >> c;
1482    while (c!=':' && !in.eof());
1483#else
1484unsigned int num;
1485do in >> c;
1486while (c!='(' && !in.eof());
1487in >> num;
1488if (num>NUMBER_DIRECTIONS) {
1489    cout << "ADOL-C error: to many directions in input\n";
1490    exit(-1);
1491}
1492do in >> c;
1493while (c!=')' && !in.eof());
1494#endif
1495    FOR_I_EQ_0_LT_NUMDIR
1496    in >> a.ADVAL_I;
1497    do in >> c;
1498    while (c!=')' && !in.eof());
1499    return in;
1500}
1501}
1502
1503/****************************************************************************/
1504#endif /* ADOLC_TAPELESS */
1505
1506/****************************************************************************/
1507/*                                                                THAT'S ALL*/
1508#endif /* __cplusplus */
1509#endif /* ADOLC_ADOUBLE_H */
Note: See TracBrowser for help on using the repository browser.