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

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

Move all external headers into a separate directory

This let's us get rid of the symlink adolc that was previously required
for building and created problems sometimes
This also adjusts all Makefiles for examples to make external builds
possible.

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

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