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

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

Merge branch 'master' of 'gitclone' into svn trunk

The following changes were merged:
commit 02942597253216bdf58b5c7f2f9a5288ce274777
Author: Max Sagebaum <sagebaum@…>
Date: Wed Apr 3 10:55:33 2013 +0200

Allow adouble objects to be initialized late

This functionality is specially required for large codes that allocate
memory using malloc/calloc and we only wish to substitute the double
datatype with adouble for tracing purposes. With malloc/calloc the
adouble constructor is not called, creating problems later in the code.
Late initialization defers some of the initialization to the first use
of an allocated space as an adouble object. This is only possible because
the class heirarchy in ADOL-C does not require polymorphism, and virtual
tables to be allocated, which can only be done in a chained superclass
constructor call.

The functionality is disabled by default and can be enabled by a configure
switch --enable-lateinit

This code is experimental and has only be tested with particular codes.

From: Max Sagebaum <sagebaum@…>
Signed-off-by: Kshitij Kulshreshtha <kshitij@…>
CC: Andrea Walther <andrea.walther@…>
CC: Nicolas Gauger <gauger@…>

commit 768851025e58a7605f0e9b6a65f51bfc3c2977f4
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Tue Apr 2 11:14:50 2013 +0200

add dvi file to gitignore

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

commit 80357742a527a4f4b9ef6323553279eb260cc7df
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Tue Apr 2 11:14:37 2013 +0200

remove as many compiler warnings from gcc -Wall as possible

only -Wunused-value remains, as we do indeed ignore some
values computed when we are in INDO, NONL_IND or
INT_FOR, INT_REV modes.

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

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