source: trunk/ADOL-C/src/adouble.h @ 71

Last change on this file since 71 was 61, checked in by awalther, 10 years ago

inclusion of error function for gcc compiler

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