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

Last change on this file since 428 was 428, checked in by kulshres, 9 years ago

get rid of fmin/fmax for doubles in traceless case

this is the traceless companion of git commit

1da1fa5bc7f5b9eada18c25026fd8ff25c8a91bb

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

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