// @HEADER
// ***********************************************************************
//
//                           Sacado Package
//                 Copyright (2007) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
// (etphipp@sandia.gov).
//
// ***********************************************************************
// @HEADER

// Extension of the RAD package (Reverse Automatic Differentiation) --
// a package specialized for function and gradient evaluations -- to
// Hessian-vector products.
// Written in 2007 by David M. Gay at Sandia National Labs, Albuquerque, NM.

#ifndef SACADO_TRAD2_H
#define SACADO_TRAD2_H

#include "Sacado_ConfigDefs.h"
#include "Sacado_trad2_Traits.hpp"

#include <stddef.h>
#include <cmath>
#include <math.h>

#if defined(RAD_DEBUG_BLOCKKEEP) && !defined(HAVE_SACADO_UNINIT)
#undef RAD_DEBUG_BLOCKKEEP
#endif

#ifdef RAD_Const_WARN	// ==> RAD_AUTO_AD_Const and RAD_DEBUG
#ifndef RAD_AUTO_AD_Const
#define RAD_AUTO_AD_Const
#endif
#ifndef RAD_DEBUG
#define RAD_DEBUG
#endif
extern "C" int RAD_Const_Warn(const void*);// outside any namespace for
					// ease in setting breakpoints
#endif // RAD_Const_WARN

#ifdef RAD_DEBUG
#include <cstdio>
#include <stdlib.h>
#endif

#ifndef RAD_AUTO_AD_Const
#ifdef RAD_DEBUG_BLOCKKEEP
#include <complex>	// must be here when SACADO_NAMESPACE is #defined
#endif
#endif

namespace Sacado {
namespace Rad2 {

// -DRAD_NO_USING_STDCC is needed, e.g., with Sun CC 5.7
#ifndef RAD_NO_USING_STDCC
  // Bring math functions into scope
  using std::exp;
  using std::log;
  using std::log10;
  using std::sqrt;
  using std::cos;
  using std::sin;
  using std::tan;
  using std::acos;
  using std::asin;
  using std::atan;
  using std::cosh;
  using std::sinh;
  using std::tanh;
  using std::abs;
  using std::fabs;
  using std::atan2;
  using std::pow;
#endif //!RAD_NO_USING_STDCC

#ifdef RAD_AUTO_AD_Const
#undef RAD_DEBUG_BLOCKKEEP
#else /*!RAD_AUTO_AD_Const*/
#ifdef RAD_DEBUG_BLOCKKEEP
#if !(RAD_DEBUG_BLOCKKEEP > 0)
#undef RAD_DEBUG_BLOCKKEEP
#else
extern "C" void _uninit_f2c(void *x, int type, long len);

template <typename T>
struct UninitType {};

template <>
struct UninitType<float> {
  static const int utype = 4;
};

template <>
struct UninitType<double> {
  static const int utype = 5;
};

template <typename T>
struct UninitType< std::complex<T> > {
  static const int utype = UninitType<T>::utype + 2;
};

#endif /*RAD_DEBUG_BLOCKKEEP > 0*/
#endif /*RAD_DEBUG_BLOCKKEEP*/
#endif /*RAD_AUTO_AD_Const*/

 class RAD_DoubleIgnore {};

 template<typename T> class
DoubleAvoid {
 public:
	typedef double	dtype;
	typedef T	ttype;
	};
 template<> class
DoubleAvoid<double> {
 public:
	typedef RAD_DoubleIgnore &dtype;
	typedef RAD_DoubleIgnore &ttype;
	};

#define Dtype typename DoubleAvoid<Double>::dtype
#define Ttype typename DoubleAvoid<Double>::ttype

 template<typename Double> class IndepADvar;
 template<typename Double> class ConstADvar;
 template<typename Double> class ConstADvari;
 template<typename Double> class ADvar;
 template<typename Double> class ADvar1;
 template<typename Double> class ADvar1g;
 template<typename Double> class ADvar1s;
 template<typename Double> class ADvar2;
 template<typename Double> class ADvar2g;
 template<typename Double> class ADvar2q;
 template<typename Double> class ADvari;
 template<typename Double> class ADvari_block;
 template<typename Double> class ADvarn;
 template<typename Double> class Derp;

 template<typename Double> struct
ADmemblock {	// We get memory in ADmemblock chunks and never give it back,
		// but reuse it once computations start anew after call(s) on
		// ADcontext::Gradcomp() or ADcontext::Weighted_Gradcomp().
	ADmemblock *next;
	Double memblk[2000];
	};

 template<typename Double> class
ADvari_block {
 public:
	typedef ADvari<Double> ADVari;
	enum { Gulp = 1021 };
	ADvari_block *next, *prev;
	ADVari **limit;
	ADVari *pADvari[Gulp];
	};

 template<typename Double> class
ADcontext {	// A singleton class: one instance in radops.c
	typedef ADmemblock<Double> ADMemblock;
	typedef ADvari <Double> ADVari;
	typedef ADvar1 <Double> ADVar1;
	typedef ADvar2 <Double> ADVar2;
	typedef ADvarn <Double> ADVarn;
	typedef ADvar1g<Double> ADVar1g;
	typedef ADvar1s<Double> ADVar1s;
	typedef ADvar2g<Double> ADVar2g;
	typedef ADvar2q<Double> ADVar2q;
	typedef ADvari_block<Double> ADVari_block;
	ADMemblock *Busy, *Free;
	char *Mbase;
	size_t Mleft;
	ADVari **Ailimit, **Ainext;
	ADVari_block *Aibusy, *Aifree;
	ADMemblock *First;
	ADVari_block *AiFirst;
	double First0[(sizeof(ADMemblock) + sizeof(double) - 1) / sizeof(double)];
	double First1[(sizeof(ADVari_block) + sizeof(double) - 1) / sizeof(double)];
	void *new_ADmemblock(size_t);
	void new_ADvari_block();
	typedef ADvar<Double> ADVar;
	typedef Derp<Double> DErp;
	int rad_need_reinit;
	size_t rad_mleft_save;
#ifdef RAD_DEBUG_BLOCKKEEP
	int rad_busy_blocks;
	ADMemblock *rad_Oldcurmb;
#endif
 public:
	static const Double One, negOne;
	ADcontext();
	void *Memalloc(size_t len);
	static void Gradcomp(int wantgrad);
	static void Hvprod(int, ADVar**, Double*, Double*);
	static void aval_reset(void);
	static void Weighted_Gradcomp(int, ADVar**, Double*);
	static inline void Gradcomp() { Gradcomp(1); }
	inline void ADvari_record(ADVari *x) {
		if (Ainext >= Ailimit)
			new_ADvari_block();
		*Ainext++ = x;
		}
	};

 template<typename Double> class
CADcontext: public ADcontext<Double> {	// for possibly constant ADvar values
 protected:
	bool fpval_implies_const;
 public:
	friend class ADvar<Double>;
	CADcontext(): ADcontext<Double>() { fpval_implies_const = false; }
	};

 template<typename Double> class
Derp {		// one derivative-propagation operation
 public:
	friend class ADvarn<Double>;
	typedef ADvari<Double> ADVari;
	static Derp *LastDerp;
	Derp *next;
	const Double *a;
	const ADVari *b;
	mutable ADVari *c;
	Derp(){};
	Derp(const ADVari *);
	Derp(const Double *, const ADVari *);
	Derp(const Double *, const ADVari *, const ADVari *);
	/* c->aval += a * b->aval; */
	};


// Now we use #define to overcome bad design in the C++ templating system

#define Ai const ADvari<Double>&
#define AI const IndepADvar<Double>&
#define T template<typename Double>
#define D Double
#define T1(f) \
T F f (AI); \
T F f (Ai);
#define T2(r,f) \
 T r f(Ai,Ai); \
 T r f(Ai,D); \
 T r f(Ai,Dtype); \
 T r f(Ai,long); \
 T r f(Ai,int); \
 T r f(D,Ai); \
 T r f(Dtype,Ai); \
 T r f(long,Ai); \
 T r f(int,Ai); \
 T r f(AI,D); \
 T r f(AI,Dtype); \
 T r f(AI,long); \
 T r f(AI,int); \
 T r f(D,AI); \
 T r f(Dtype,AI); \
 T r f(long,AI); \
 T r f(int,AI); \
 T r f(Ai,AI);\
 T r f(AI,Ai);\
 T r f(AI,AI);

#define F ADvari<Double>&
T2(F, operator+)
T2(F, operator-)
T2(F, operator*)
T2(F, operator/)
T2(F, atan2)
T2(F, pow)
T2(F, max)
T2(F, min)
T2(int, operator<)
T2(int, operator<=)
T2(int, operator==)
T2(int, operator!=)
T2(int, operator>=)
T2(int, operator>)
T1(operator+)
T1(operator-)
T1(abs)
T1(acos)
T1(acosh)
T1(asin)
T1(asinh)
T1(atan)
T1(atanh)
T1(cos)
T1(cosh)
T1(exp)
T1(log)
T1(log10)
T1(sin)
T1(sinh)
T1(sqrt)
T1(tan)
T1(tanh)
T1(fabs)

T F copy(AI);
T F copy(Ai);

#undef F
#undef T2
#undef T1
#undef D
#undef T
#undef AI
#undef Ai

} /* namespace Rad2 */
} /* namespace Sacado */
#define SNS Sacado::Rad2
namespace std {	// Moved here from bottom for use in testing nesting of Rad with itself.
  using SNS::exp;
  using SNS::log;
  using SNS::log10;
  using SNS::sqrt;
  using SNS::cos;
  using SNS::sin;
  using SNS::tan;
  using SNS::acos;
  using SNS::asin;
  using SNS::atan;
  using SNS::cosh;
  using SNS::sinh;
  using SNS::tanh;
  using SNS::abs;
  using SNS::fabs;
  using SNS::atan2;
  using SNS::pow;
}
#undef SNS
namespace Sacado {
namespace Rad2 {

 template<typename Double>ADvari<Double>& ADf1(Double f, Double g, const IndepADvar<Double> &x);
 template<typename Double>ADvari<Double>& ADf2(Double f, Double gx, Double gy,
	const IndepADvar<Double> &x, const IndepADvar<Double> &y);
 template<typename Double>ADvari<Double>& ADfn(Double f, int n,
	const IndepADvar<Double> *x, const Double *g);

 template<typename Double> IndepADvar<Double>& ADvar_operatoreq(IndepADvar<Double>*,
	const ADvari<Double>&);
 template<typename Double> ADvar<Double>& ADvar_operatoreq(ADvar<Double>*, const ADvari<Double>&);
 template<typename Double> void AD_Const(const IndepADvar<Double>&);
 template<typename Double> void AD_Const1(Double*, const IndepADvar<Double>&);
 template<typename Double> ADvari<Double>& ADf1(Double, Double, const ADvari<Double>&);
 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
	const ADvari<Double>&, const ADvari<Double>&);
 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
	const IndepADvar<Double>&, const ADvari<Double>&);
 template<typename Double> ADvari<Double>& ADf2(Double, Double, Double,
	const ADvari<Double>&, const IndepADvar<Double>&);
 template<typename Double> Double val(const ADvari<Double>&);

 enum Advari_Opclass {
	Hv_const,
	Hv_copy,
	Hv_binary,
	Hv_unary,
	Hv_negate,
	Hv_plusLR,
	Hv_minusLR,
	Hv_timesL,
	Hv_timesLR,
	Hv_quotLR,
	Hv_nary
	};

 template<typename Double> ADvari<Double>&
ADf1(Double f, Double g, Double h, const ADvari<Double> &x);

 template<typename Double> ADvari<Double>&
ADf2(Double f, Double gx, Double gy, Double hxx,
			Double hxy, Double hyy, const ADvari<Double> &x, const ADvari<Double> &y);

 template<typename Double> ADvari<Double>&
ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h);

 template<typename Double> class
ADvari {	// implementation of an ADvar
 public:
	typedef Double value_type;
	typedef IndepADvar<Double> IndepADVar;
#ifdef RAD_AUTO_AD_Const
	friend class IndepADvar<Double>;
#ifdef RAD_Const_WARN
	mutable const IndepADVar *padv;
#else
 protected:
	mutable const IndepADVar *padv;
#endif //RAD_Const_WARN
 private:
	ADvari *Next;
	static ADvari *First_ADvari, **Last_ADvari;
 public:
	inline void ADvari_padv(const IndepADVar *v) {
		*Last_ADvari = this;
		Last_ADvari = &this->Next;
		this->padv = v;
		}
#endif //RAD_AUTO_AD_Const
#ifdef RAD_DEBUG
	int gcgen;
	int opno;
	static int gcgen_cur, last_opno, zap_gcgen, zap_gcgen1, zap_opno;
	static FILE *debug_file;
#endif
	static ADcontext<Double> adc;
	Advari_Opclass opclass;
	Double Val;		// result of this operation
	mutable Double aval;	// adjoint -- partial of final result w.r.t. this Val
	mutable Double dO;	// deriv of op w.r.t. t in x + t*p
	mutable Double aO;	// adjoint (in Hv computation) of op
	mutable Double adO;	// adjoint (in Hv computation) of dO
	void *operator new(size_t len) {
#ifdef RAD_DEBUG
		ADvari *rv = (ADvari*)ADvari::adc.Memalloc(len);
		rv->gcgen = gcgen_cur;
		rv->opno = ++last_opno;
		if (last_opno == zap_opno && gcgen_cur == zap_gcgen)
			printf("");
		return rv;
#else
		return ADvari::adc.Memalloc(len);
#endif
		}
	void operator delete(void*) {} /*Should never be called.*/
	inline ADvari(Advari_Opclass oc, Double t):
		opclass(oc), Val(t), aval(0.), dO(0.)
		{ if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
	inline ADvari(Advari_Opclass oc, Double t, Double ta):
		opclass(oc), Val(t), aval(ta), dO(0.)
		{ if (oc != Hv_const) ADvari::adc.ADvari_record(this); }
 private:
	inline ADvari(): Val(0.), aval(0.), dO(0.) {}	// prevent construction without value (?)
 public:
	friend class ConstADvari<Double>;
#ifdef RAD_AUTO_AD_Const
	friend class ADcontext<Double>;
	friend class ADvar<Double>;
	friend class ADvar1<Double>;
	friend class ADvar1s<Double>;
	friend class ADvar2<Double>;
	friend class ADvar2q<Double>;
	friend class ConstADvar<Double>;
	ADvari(const IndepADVar *, Double);
	ADvari(const IndepADVar *, Double, Double);
	ADvari(const IndepADVar *, Double, Double, int);
#endif
#define F friend
#define R ADvari&
#define Ai const ADvari&
#define T1(r,f) F r f <>(Ai);
#define T2(r,f) \
F r f <>(Ai,Ai); \
F r f <>(Ttype,Ai); \
F r f <>(Ai,Ttype); \
F r f <>(double,Ai); \
F r f <>(Ai,double); \
F r f <>(long,Ai); \
F r f <>(Ai,long); \
F r f <>(int,Ai); \
F r f <>(Ai,int);
	T1(R,operator+)
	T2(R,operator+)
	T1(R,operator-)
	T2(R,operator-)
	T2(R,operator*)
	T2(R,operator/)
	T1(R,abs)
	T1(R,acos)
	T1(R,acosh)
	T1(R,asin)
	T1(R,asinh)
	T1(R,atan)
	T1(R,atanh)
	T2(R,atan2)
	T2(R,max)
	T2(R,min)
	T1(R,cos)
	T1(R,cosh)
	T1(R,exp)
	T1(R,log)
	T1(R,log10)
	T2(R,pow)
	T1(R,sin)
	T1(R,sinh)
	T1(R,sqrt)
	T1(R,tan)
	T1(R,tanh)
	T1(R,fabs)
	T1(R,copy)
	T2(int,operator<)
	T2(int,operator<=)
	T2(int,operator==)
	T2(int,operator!=)
	T2(int,operator>=)
	T2(int,operator>)
#undef T2
#undef T1
#undef Ai
#undef R
#undef F

	friend ADvari& ADf1<>(Double f, Double g, Double h, const ADvari &x);
	friend ADvari& ADf2<>(Double f, Double gx, Double gy, Double hxx,
			Double hxy, Double hyy, const ADvari &x, const ADvari &y);
	friend ADvari& ADfn<>(Double f, int n, const IndepADVar *x,
			const Double *g, const Double *h);

	inline operator Double() { return this->Val; }
	inline operator Double() const { return this->Val; }
	};

 template<typename Double> class
ADvar1: public ADvari<Double> {	// simplest unary ops
 public:
	typedef ADvari<Double> ADVari;
	Derp<Double> d;
	ADvar1(Advari_Opclass oc, Double val1, const Double *a1, const ADVari *c1):
		ADVari(oc,val1), d(a1,this,c1) {}
#ifdef RAD_AUTO_AD_Const
	typedef typename ADVari::IndepADVar IndepADVar;
	typedef ADvar<Double> ADVar;
	ADvar1(const IndepADVar*, const IndepADVar&);
	ADvar1(const IndepADVar*, const ADVari&);
	ADvar1(Advari_Opclass oc, const Double val1, const Double *a1,
		const ADVari *c1, const ADVar *v):
			ADVari(oc,val1), d(a1,this,c1) {
		c1->padv = 0;
		this->ADvari_padv(v);
		}
#endif
	};


 template<typename Double> class
ConstADvari: public ADvari<Double> {
 private:
	ConstADvari *prevcad;
	ConstADvari() {};	// prevent construction without value (?)
	static ConstADvari *lastcad;
 public:
	typedef ADvari<Double> ADVari;
	typedef Derp<Double> DErp;
	static CADcontext<Double> cadc;
	inline void *operator new(size_t len) { return ConstADvari::cadc.Memalloc(len); }
	inline ConstADvari(Double t): ADVari(Hv_copy, t) { prevcad = lastcad; lastcad = this; }
	static void aval_reset(void);
	};


 template<typename Double> class
IndepADvar {		// an independent ADvar
 private:
	IndepADvar& operator=(IndepADvar&x) {
		/* private to prevent assignment */
#ifdef RAD_AUTO_AD_Const
		if (cv)
			cv->padv = 0;
		return new ADvar1<Double>(this,x);
#elif defined(RAD_EQ_ALIAS)
		this->cv = x.cv;
		return *this;
#else
		return ADvar_operatoreq(this,*x.cv);
#endif //RAD_AUTO_AD_Const
		}
 protected:
	static void AD_Const(const IndepADvar&);
	mutable ADvari<Double> *cv;
 public:
	typedef Double value_type;
	friend class ADvar<Double>;
	friend class ADcontext<Double>;
	friend class ADvar1<Double>;
	friend class ADvarn<Double>;
	typedef ADvari<Double> ADVari;
	typedef ADvar<Double> ADVar;
	IndepADvar(Ttype);
	IndepADvar(double);
	IndepADvar(int);
	IndepADvar(long);
	IndepADvar& operator= (Double);
#ifdef RAD_AUTO_AD_Const
	friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
	inline IndepADvar(const IndepADvar &x) { cv = x.cv ? new ADvar1<Double>(this, x) : 0; };
	inline IndepADvar(ADVari *ncv) { cv = ncv; }
	inline IndepADvar() { cv = 0; }
	inline ~IndepADvar() {
					if (cv)
						cv->padv = 0;
					}
#else
	inline IndepADvar() {
#ifndef RAD_EQ_ALIAS
		cv = 0;
#endif
		}
	inline ~IndepADvar() {}
	friend IndepADvar& ADvar_operatoreq<>(IndepADvar*, const ADVari&);
#endif

#ifdef RAD_Const_WARN
	inline operator ADVari&() const {
		ADVari *tcv = this->cv;
		if (tcv->opno < 0)
			RAD_Const_Warn(tcv);
		return *tcv;
		}
	inline operator ADVari*() const {
		ADVari *tcv = this->cv;
		if (tcv->opno < 0)
			RAD_Const_Warn(tcv);
		return tcv;
		}
#else //RAD_Const_WARN
	inline operator ADVari&() const { return *this->cv; }
	inline operator ADVari*() const { return this->cv; }
#endif //RAD_Const_WARN

	Double val() const { return cv->Val; }
	Double adj() const { return cv->aval; }

	friend void AD_Const1<>(Double*, const IndepADvar&);

	friend ADVari& ADf1<>(Double, Double, const IndepADvar&);
	friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const IndepADvar&);
	friend ADVari& ADf2<>(Double, Double, Double, const ADVari&, const IndepADvar&);
	friend ADVari& ADf2<>(Double, Double, Double, const IndepADvar&, const ADVari&);

	static inline void Gradcomp(int wantgrad)
				{ ADcontext<Double>::Gradcomp(wantgrad); }
	static inline void Gradcomp()
				{ ADcontext<Double>::Gradcomp(1); }
	static inline void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
				{ ADcontext<Double>::Hvprod(n, vp, v, hv); }
	static inline void aval_reset() { ConstADvari<Double>::aval_reset(); }
	static inline void Weighted_Gradcomp(int n, ADVar **v, Double *w)
				{ ADcontext<Double>::Weighted_Gradcomp(n, v, w); }

	/* We use #define to deal with bizarre templating rules that apparently */
	/* require us to spell the some conversion explicitly */


#define Ai const ADVari&
#define AI const IndepADvar&
#define D Double
#define T2(r,f) \
 r f <>(AI,AI);\
 r f <>(Ai,AI);\
 r f <>(AI,Ai);\
 r f <>(Ttype,AI);\
 r f <>(double,AI);\
 r f <>(long,AI);\
 r f <>(int,AI);\
 r f <>(AI,Ttype);\
 r f <>(AI,double);\
 r f <>(AI,long);\
 r f <>(AI,int);
#define T1(f) friend ADVari& f<> (AI);

#define F friend ADVari&
T2(F, operator+)
T2(F, operator-)
T2(F, operator*)
T2(F, operator/)
T2(F, atan2)
T2(F, max)
T2(F, min)
T2(F, pow)
#undef F
#define F friend int
T2(F, operator<)
T2(F, operator<=)
T2(F, operator==)
T2(F, operator!=)
T2(F, operator>=)
T2(F, operator>)

T1(operator+)
T1(operator-)
T1(abs)
T1(acos)
T1(acosh)
T1(asin)
T1(asinh)
T1(atan)
T1(atanh)
T1(cos)
T1(cosh)
T1(exp)
T1(log)
T1(log10)
T1(sin)
T1(sinh)
T1(sqrt)
T1(tan)
T1(tanh)
T1(fabs)
T1(copy)

#undef F
#undef T1
#undef T2
#undef D
#undef AI
#undef Ai

	};

 template<typename Double> class
ADvar: public IndepADvar<Double> {	// an "active" variable
 public:
        //! Turn ADvar into a meta-function class usable with mpl::apply
        template <typename U> struct apply { typedef ADvar<U> type; };

	typedef IndepADvar<Double> IndepADVar;
	typedef ADvari<Double> ADVari;
	typedef ConstADvari<Double> ConstADVari;
 private:
	void ADvar_ctr(Double d) {
		ADVari *x;
		if (ConstADVari::cadc.fpval_implies_const)
			x = new ConstADVari(d);
		else {
#ifdef RAD_AUTO_AD_Const
			x = new ADVari((IndepADVar*)this, d);
			x->ADvari_padv(this);
#else
			x = new ADVari(Hv_const, d);
#endif
			}
		this->cv = x;
		}
 public:
	friend class ADvar1<Double>;
	typedef ADvar1<Double> ADVar1;
	ADvar() { /* cv = 0; */ }
	ADvar(Ttype d)  { ADvar_ctr(d); }
	ADvar(double i) { ADvar_ctr(Double(i)); }
	ADvar(int i)	{ ADvar_ctr(Double(i)); }
	ADvar(long i)	{ ADvar_ctr(Double(i)); }
	inline ~ADvar() {}
#ifdef RAD_AUTO_AD_Const
	ADvar(IndepADVar &x) {
		this->cv = x.cv ? new ADVar1(this, x) : 0;
		}
	ADvar(ADVari &x) :IndepADVar(&x) { x.ADvari_padv((IndepADVar*)this);}
	inline ADvar& operator=(IndepADVar &x) {
		if (this->cv)
			this->cv->padv = 0;
		this->cv = new ADVar1(this,x);
		return *this;
		}
	inline ADvar& operator=(ADVari &x) {
		if (this->cv)
			this->cv->padv = 0;
		this->cv = new ADVar1(this, x);
		return *this;
		}
#else /*!RAD_AUTO_AD_Const*/
	friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
#ifdef RAD_EQ_ALIAS
	/* allow aliasing v and w after "v = w;" */
	inline ADvar(const IndepADVar &x) { this->cv = (ADVari*)x.cv; }
	inline ADvar(const ADVari &x) { this->cv = (ADVari*)&x; }
	inline ADvar& operator=(IndepADVar &x) { this->cv = (ADVari*)x.cv; return *this; }
	inline ADvar& operator=(const ADVari &x) { this->cv = (ADVari*)&x; return *this; }
#else /*!RAD_EQ_ALIAS*/
	friend ADvar& ADvar_operatoreq<>(ADvar*, const ADVari&);
	ADvar(const IndepADVar &x) { this->cv = x.cv ? new ADVar1(x.cv->Val, &this->cv->adc.One, x.cv) : 0; }
	ADvar(const ADvar&x) { this->cv = x.cv ?
		new ADVar1(Hv_copy, x.cv->Val, &this->cv->adc.One, (ADVari*)x.cv) : 0; }
	ADvar(const ADVari &x) { this->cv = new ADVar1(Hv_copy, x.Val, &this->cv->adc.One, &x); }
	inline ADvar& operator=(IndepADVar &x) { return ADvar_operatoreq(this,*x.cv); };
	inline ADvar& operator=(const ADVari &x) { return ADvar_operatoreq(this,x); };
#endif /* RAD_EQ_ALIAS */
#endif /* RAD_AUTO_AD_Const */
	ADvar& operator=(Double);
	ADvar& operator+=(const ADVari&);
	ADvar& operator+=(Double);
	ADvar& operator-=(const ADVari&);
	ADvar& operator-=(Double);
	ADvar& operator*=(const ADVari&);
	ADvar& operator*=(Double);
	ADvar& operator/=(const ADVari&);
	ADvar& operator/=(Double);
	inline static bool get_fpval_implies_const(void)
		{ return ConstADVari::cadc.fpval_implies_const; }
	inline static void set_fpval_implies_const(bool newval)
		{ ConstADVari::cadc.fpval_implies_const = newval; }
	inline static bool setget_fpval_implies_const(bool newval) {
		bool oldval = ConstADVari::cadc.fpval_implies_const;
		ConstADVari::cadc.fpval_implies_const = newval;
		return oldval;
		}
	static inline void Gradcomp(int wantgrad)
				{ ADcontext<Double>::Gradcomp(wantgrad); }
	static inline void Gradcomp()
				{ ADcontext<Double>::Gradcomp(1); }
	static inline void aval_reset() { ConstADVari::aval_reset(); }
	static inline void Weighted_Gradcomp(int n, ADvar **v, Double *w)
				{ ADcontext<Double>::Weighted_Gradcomp(n, v, w); }
	};

template<typename Double>
 inline void AD_Const1(Double *notused, const IndepADvar<Double>&v)
{ IndepADvar<Double>::AD_Const(v); }

template<typename Double>
 inline void AD_Const(const IndepADvar<Double>&v) { AD_Const1((Double*)0, v); }

 template<typename Double> class
ConstADvar: public ADvar<Double> {
 public:
	typedef ADvar<Double> ADVar;
	typedef typename ADVar::ADVari ADVari;
	typedef typename ADVar::ConstADVari ConstADVari;
	typedef Derp<Double> DErp;
	typedef typename ADVar::IndepADVar IndepADVar;
 private: // disable op=
	ConstADvar& operator+=(ADVari&);
	ConstADvar& operator+=(Double);
	ConstADvar& operator-=(ADVari&);
	ConstADvar& operator-=(Double);
	ConstADvar& operator*=(ADVari&);
	ConstADvar& operator*=(Double);
	ConstADvar& operator/=(ADVari&);
	ConstADvar& operator/=(Double);
	void ConstADvar_ctr(Double);
 public:
	ConstADvar(Ttype d)	{ ConstADvar_ctr(d); }
	ConstADvar(double i)	{ ConstADvar_ctr(Double(i)); }
	ConstADvar(int i)	{ ConstADvar_ctr(Double(i)); }
	ConstADvar(long i)	{ ConstADvar_ctr(Double(i)); }
	ConstADvar(const IndepADVar&);
	ConstADvar(const ConstADvar&);
	ConstADvar(const ADVari&);
	inline ~ConstADvar() {}
#ifdef RAD_NO_CONST_UPDATE
 private:
#endif
	ConstADvar();
	inline ConstADvar& operator=(Double d) { this->cv->Val = d; return *this; }
	inline ConstADvar& operator=(ADVari& d) { this->cv->Val = d.Val; return *this; }
 };

 template<typename Double> class
ADvar1s: public ADvar1<Double> { // unary ops with partials
 public:
	typedef ADvar1<Double> ADVar1;
	typedef typename ADVar1::ADVari ADVari;
	Double pL;	// deriv of op w.r.t. left operand L
	ADvar1s(Double val1, Double a1, const ADVari *c1):
		ADVar1(Hv_timesL,val1,&pL,c1), pL(a1) {}
#ifdef RAD_AUTO_AD_Const
	Double a;
	typedef typename ADVar1::ADVar ADVar;
	ADvar1s(Double val1, Double a1, const ADVari *c1, const ADVar *v):
		ADVar1(Hv_timesL,val1,&a,c1,v), a(a1) {}
#endif
	};

 template<typename Double> class
ADvar1g: public ADvar1<Double> { // unary ops with partials
 public:
	typedef ADvar1<Double> ADVar1;
	typedef typename ADVar1::ADVari ADVari;
	Double pL;	// deriv of op w.r.t. left operand L
	Double pL2;	// partial of op w.r.t. L,L
	ADvar1g(Double val1, Double d1, Double d2, const ADVari *c1):
		ADVar1(Hv_unary,val1,&pL,c1), pL(d1), pL2(d2) {}
	};

 template<typename Double> class
ADvar2: public ADvari<Double> {	// basic binary ops
 public:
	typedef ADvari<Double> ADVari;
	typedef Derp<Double> DErp;
	DErp dL, dR;
	ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
		const ADVari *Rcv, const Double *Rc):
			ADVari(oc,val1) {
		dR.next = DErp::LastDerp;
		dL.next = &dR;
		DErp::LastDerp = &dL;
		dL.a = Lc;
		dL.c = (ADVari*)Lcv;
		dR.a = Rc;
		dR.c = (ADVari*)Rcv;
		dL.b = dR.b = this;
		}
#ifdef RAD_AUTO_AD_Const
	typedef ADvar<Double> ADVar;
	ADvar2(Advari_Opclass oc, Double val1, const ADVari *Lcv, const Double *Lc,
		const ADVari *Rcv, const Double *Rc, ADVar *v):
			ADVari(oc,val1) {
		dR.next = DErp::LastDerp;
		dL.next = &dR;
		DErp::LastDerp = &dL;
		dL.a = Lc;
		dL.c = Lcv;
		dR.a = Rc;
		dR.c = Rcv;
		dL.b = dR.b = this;
		Lcv->padv = 0;
		this->ADvari_padv(v);
		}
#endif
	};

 template<typename Double> class
ADvar2q: public ADvar2<Double> { // binary ops with partials
 public:
	typedef ADvar2<Double> ADVar2;
	typedef typename ADVar2::ADVari ADVari;
	typedef typename ADVar2::DErp DErp;
	Double pL;	// deriv of op w.r.t. left operand L
	Double pR;	// deriv of op w.r.t. right operand R
	Double pLR;	// second partial w.r.t. L,R
	Double pR2;	// second partial w.r.t. R,R
	ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
		const ADVari *Lcv, const ADVari *Rcv):
			ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR),
			pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
#ifdef RAD_AUTO_AD_Const
	typedef typename ADVar2::ADVar ADVar;
	ADvar2q(Double val1, Double Lp, Double Rp, Double LR, Double R2,
		const ADVari *Lcv, const ADVari *Rcv, const ADVar *v):
			ADVar2(Hv_quotLR,val1,Lcv,&pL,Rcv,&pR,v),
			pL(Lp), pR(Rp), pLR(LR), pR2(R2) {}
#endif
	};

 template<typename Double> class
ADvar2g: public ADvar2<Double> { // general binary ops with partials
 public:
	typedef ADvar2<Double> ADVar2;
	typedef typename ADVar2::ADVari ADVari;
	Double pL;	// deriv of op w.r.t. left operand L
	Double pR;	// deriv of op w.r.t. right operand R
	Double pL2;	// second partial w.r.t. L,L
	Double pLR;	// second partial w.r.t. L,R
	Double pR2;	// second partial w.r.t. R,R
	ADvar2g(Double val1, Double Lp, Double Rp, Double L2, Double LR, Double R2,
		const ADVari *Lcv, const ADVari *Rcv):
			ADVar2(Hv_binary,val1,Lcv,&pL,Rcv,&pR),
			pL(Lp), pR(Rp), pL2(L2), pLR(LR), pR2(R2) { }
	};

 template<typename Double> class
ADvarn: public ADvari<Double> { // n-ary ops with partials g and
				// 2nd partials h (lower triangle, rowwise)
 public:
	typedef ADvari<Double> ADVari;
	typedef ADvar<Double> ADVar;
	typedef typename ADVari::IndepADVar IndepADVar;
	typedef Derp<Double> DErp;
	int n;
	Double *G, *H;
	DErp *D;
	ADvarn(Double val1, int n1, const IndepADVar *x, const Double *g, const Double *h):
			ADVari(Hv_nary,val1), n(n1) {
		DErp *d1, *dlast;
		Double *a1;
		int i, nh;

		a1 = G = (Double*)ADVari::adc.Memalloc(n1*sizeof(*g));
		d1 = D = (DErp*)ADVari::adc.Memalloc(n1*sizeof(DErp));
		dlast = DErp::LastDerp;
		for(i = 0; i < n1; i++, d1++) {
			d1->next = dlast;
			dlast = d1;
			a1[i] = g[i];
			d1->a = &a1[i];
			d1->b = this;
			d1->c = x[i].cv;
			}
		DErp::LastDerp = dlast;
		nh = (n1*(n1+1)) >> 1;
		a1 = H = (double*)ADVari::adc.Memalloc(nh * sizeof(*h));
		for(i = 0; i < nh; i++)
			a1[i] = h[i];
		}
	};

template<typename Double>
 inline ADvari<Double>& operator+(const ADvari<Double> &T) { return *(ADvari<Double>*)&T; }

template<typename Double>
 inline int operator<(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val < R.Val; }
template<typename Double>
 inline int operator<(const ADvari<Double> &L, Double R) { return L.Val < R; }
template<typename Double>
 inline int operator<(Double L, const ADvari<Double> &R) { return L < R.Val; }

template<typename Double>
 inline int operator<=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val <= R.Val; }
template<typename Double>
 inline int operator<=(const ADvari<Double> &L, Double R) { return L.Val <= R; }
template<typename Double>
 inline int operator<=(Double L, const ADvari<Double> &R) { return L <= R.Val; }

template<typename Double>
 inline int operator==(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val == R.Val; }
template<typename Double>
 inline int operator==(const ADvari<Double> &L, Double R) { return L.Val == R; }
template<typename Double>
 inline int operator==(Double L, const ADvari<Double> &R) { return L == R.Val; }

template<typename Double>
 inline int operator!=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val != R.Val; }
template<typename Double>
 inline int operator!=(const ADvari<Double> &L, Double R) { return L.Val != R; }
template<typename Double>
 inline int operator!=(Double L, const ADvari<Double> &R) { return L != R.Val; }

template<typename Double>
 inline int operator>=(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val >= R.Val; }
template<typename Double>
 inline int operator>=(const ADvari<Double> &L, Double R) { return L.Val >= R; }
template<typename Double>
 inline int operator>=(Double L, const ADvari<Double> &R) { return L >= R.Val; }

template<typename Double>
 inline int operator>(const ADvari<Double> &L, const ADvari<Double> &R) { return L.Val > R.Val; }
template<typename Double>
 inline int operator>(const ADvari<Double> &L, Double R) { return L.Val > R; }
template<typename Double>
 inline int operator>(Double L, const ADvari<Double> &R) { return L > R.Val; }

template<typename Double>
 inline void *ADcontext<Double>::Memalloc(size_t len) {
		if (Mleft >= len)
			return Mbase + (Mleft -= len);
		return new_ADmemblock(len);
		}

template<typename Double>
 inline Derp<Double>::Derp(const ADVari *c1): c((ADVari*)c1) {
		next = LastDerp;
		LastDerp = this;
		}

template<typename Double>
 inline Derp<Double>::Derp(const Double *a1, const ADVari *c1): a(a1), c((ADVari*)c1) {
		next = LastDerp;
		LastDerp = this;
		}

template<typename Double>
 inline Derp<Double>::Derp(const Double *a1, const ADVari *b1, const ADVari *c1):
	a(a1), b(b1), c((ADVari*)c1) {
		next = LastDerp;
		LastDerp = this;
		}

/**** radops ****/

template<typename Double> Derp<Double> *Derp<Double>::LastDerp = 0;
template<typename Double> ADcontext<Double> ADvari<Double>::adc;
template<typename Double> const Double ADcontext<Double>::One = 1.;
template<typename Double> const Double ADcontext<Double>::negOne = -1.;
template<typename Double> CADcontext<Double> ConstADvari<Double>::cadc;
template<typename Double> ConstADvari<Double> *ConstADvari<Double>::lastcad;

#ifdef RAD_AUTO_AD_Const
template<typename Double> ADvari<Double>*   ADvari<Double>::First_ADvari;
template<typename Double> ADvari<Double>**  ADvari<Double>::Last_ADvari = &ADvari<Double>::First_ADvari;
#endif

#ifdef RAD_DEBUG
#ifndef RAD_DEBUG_gcgen1
#define RAD_DEBUG_gcgen1 -1
#endif
template<typename Double> int ADvari<Double>::gcgen_cur;
template<typename Double> int ADvari<Double>::last_opno;
template<typename Double> int ADvari<Double>::zap_gcgen;
template<typename Double> int ADvari<Double>::zap_gcgen1 = RAD_DEBUG_gcgen1;
template<typename Double> int ADvari<Double>::zap_opno;
template<typename Double> FILE *ADvari<Double>::debug_file;
#endif


 template<typename Double>
ADcontext<Double>::ADcontext()
{
	ADVari_block *fb;

	First = (ADMemblock*)First0;
	First->next = 0;
	Busy = First;
	Free = 0;
	Mbase = (char*)First->memblk;
	Mleft = sizeof(First->memblk);
	AiFirst = Aibusy = fb = (ADVari_block*)First1;
	Aifree = 0;
	Ainext = fb->pADvari;
	fb->next = fb->prev = 0;
	fb->limit = Ailimit = fb->pADvari + ADVari_block::Gulp;
	rad_need_reinit = 0;
#ifdef RAD_DEBUG_BLOCKKEEP
	rad_busy_blocks = 0;
	rad_mleft_save = 0;
	rad_Oldcurmb = 0;
#endif
	}

template<typename Double> void*
ADcontext<Double>::new_ADmemblock(size_t len)
{
	ADMemblock *mb, *mb0, *mb1, *mbf, *x;
	ADVari_block *b;
#ifdef RAD_AUTO_AD_Const
	ADVari *a, *anext;
	IndepADvar<Double> *v;
#ifdef RAD_Const_WARN
	ADVari *cv;
	int i, j;
#endif
#endif /*RAD_AUTO_AD_Const*/

	if ((rad_need_reinit & 1) && this == &ADVari::adc) {
		rad_need_reinit &= ~1;
		DErp::LastDerp = 0;
		Aibusy = b = AiFirst;
		Aifree = b->next;
		b->next = b->prev = 0;
		Ailimit = b->limit = (Ainext = b->pADvari) + ADVari_block::Gulp;
#ifdef RAD_DEBUG_BLOCKKEEP
		Mleft = rad_mleft_save;
		if (Mleft < sizeof(First->memblk))
			_uninit_f2c(Mbase + Mleft,
				UninitType<Double>::utype,
			 	(sizeof(First->memblk) - Mleft)
				/sizeof(typename Sacado::ValueType<Double>::type));
		if ((mb = Busy->next)) {
			if (!(mb0 = rad_Oldcurmb))
				mb0 = (ADMemblock*)First0;
			for(;; mb = mb->next) {
				_uninit_f2c(mb->memblk,
					UninitType<Double>::utype,
					sizeof(First->memblk)
					/sizeof(typename Sacado::ValueType<Double>::type));
				if (mb == mb0)
					break;
				}
			}
		rad_Oldcurmb = Busy;
		if (rad_busy_blocks >= RAD_DEBUG_BLOCKKEEP) {
			rad_busy_blocks = 0;
			rad_Oldcurmb = 0;
			mb0 = (ADMemblock*)First0;
			mbf =  Free;
			for(mb = Busy; mb != mb0; mb = mb1) {
				mb1 = mb->next;
				mb->next = mbf;
				mbf = mb;
				}
			Free = mbf;
			Busy = mb;
			Mbase = (char*)First->memblk;
			Mleft = sizeof(First->memblk);
			}

#else /* !RAD_DEBUG_BLOCKKEEP */

		mb0 = First;
		mbf =  Free;
		for(mb = Busy; mb != mb0; mb = mb1) {
			mb1 = mb->next;
			mb->next = mbf;
			mbf = mb;
			}
		Free = mbf;
		Busy = mb;
		Mbase = (char*)First->memblk;
		Mleft = sizeof(First->memblk);
#ifdef RAD_AUTO_AD_Const
		if (ADVari::adc.rad_need_reinit & 2) {
			ADVari::adc.rad_need_reinit &= ~2;
			*ADVari::Last_ADvari = 0;
			ADVari::Last_ADvari = &ADVari::First_ADvari;
			if ((anext = ADVari::First_ADvari)) {
				while((a = anext)) {
					anext = a->Next;
					if ((v = (IndepADvar<Double> *)a->padv)) {
#ifdef RAD_Const_WARN
						if ((i = a->opno) > 0)
							i = -i;
						j = a->gcgen;
						v->cv = cv = new ADVari(v, a->Val, a->aval);
						cv->opno = i;
						cv->gcgen = j;
#else
						v->cv = new ADVari(v, a->Val, a->aval);
#endif
						}
					}
				}
			}
#endif /*RAD_AUTO_AD_Const*/
#endif /*RAD_DEBUG_BLOCKKEEP*/
		if (Mleft >= len)
			return Mbase + (Mleft -= len);
		}

	if ((x = Free))
		Free = x->next;
	else
		x = new ADMemblock;
#ifdef RAD_DEBUG_BLOCKKEEP
	rad_busy_blocks++;
#endif
	x->next = Busy;
	Busy = x;
	return (Mbase = (char*)x->memblk) +
		(Mleft = sizeof(First->memblk) - len);
	}

template<typename Double> void
ADcontext<Double>::new_ADvari_block()
{
	ADVari_block *ob, *nb;
	ob = Aibusy;
	ob->limit = Ailimit;	// should be unnecessary, but harmless
	if ((nb = Aifree))
		Aifree = nb->next;
	else
		nb = new ADVari_block;
	Aibusy = Aibusy->next = nb;
	nb->limit = Ailimit = (Ainext = nb->pADvari) + ADVari_block::Gulp;
	ob->next = nb;
	nb->prev = ob;
	nb->next = 0;
	}

template<typename Double> void
ADcontext<Double>::Gradcomp(int wantgrad)
{
	DErp *d;

	if (ADVari::adc.rad_need_reinit) {
		for(d = DErp::LastDerp; d; d = d->next)
			d->c->aval = 0;
		}
	if (!(ADVari::adc.rad_need_reinit & 1)) {
		ADVari::adc.rad_need_reinit = 1;
		ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
		ADVari::adc.Mleft = 0;
		}
#ifdef RAD_DEBUG
	if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
		const char *fname;
		if (!(fname = getenv("RAD_DEBUG_FILE")))
			fname = "rad_debug.out";
		else if (!*fname)
			fname = 0;
		if (fname)
			ADVari::debug_file = fopen(fname, "w");
		ADVari::zap_gcgen1 = -1;
		}
#endif
	if ((d = DErp::LastDerp) != 0) {
#ifdef RAD_AUTO_AD_Const
		 ADVari::adc.rad_need_reinit |= 2;
#endif /*RAD_AUTO_AD_Const*/
		if (wantgrad) {
			d->b->aval = 1;
#ifdef RAD_DEBUG
			if (ADVari::debug_file)
				do {
					fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
						d->c->opno, d->b->opno, d->c->aval,
						*d->a, d->b->aval);
					d->c->aval += *d->a * d->b->aval;
					fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
					fflush(ADVari::debug_file);
					} while((d = d->next));
			else
#endif
			do d->c->aval += *d->a * d->b->aval;
			while((d = d->next));
			DErp::LastDerp = 0;
			}
		}
#ifdef RAD_DEBUG
	if (ADVari::debug_file) {
		fclose(ADVari::debug_file);
		ADVari::debug_file = 0;
		}
#endif //RAD_DEBUG
#ifdef RAD_DEBUG
	ADVari::gcgen_cur++;
	ADVari::last_opno = 0;
#endif
	}

template<typename Double> void
ADcontext<Double>::Weighted_Gradcomp(int n, ADVar **V, Double *w)
{
	DErp *d;
	int i;
#ifdef RAD_Const_WARN
	ADVari *cv;
	int j;
#endif
#ifdef RAD_AUTO_AD_Const
	ADVari *a, *anext;
	IndepADvar<Double> *v;
#endif /*RAD_AUTO_AD_Const*/

	if (ADVari::adc.rad_need_reinit) {
		for(d = DErp::LastDerp; d; d = d->next)
			d->c->aval = 0;
		}
	if (!(ADVari::adc.rad_need_reinit & 1)) {
		ADVari::adc.rad_need_reinit = 1;
		ADVari::adc.rad_mleft_save = ADVari::adc.Mleft;
		ADVari::adc.Mleft = 0;
		}
#ifdef RAD_DEBUG
	if (ADVari::gcgen_cur == ADVari::zap_gcgen1) {
		const char *fname;
		if (!(fname = getenv("RAD_DEBUG_FILE")))
			fname = "rad_debug.out";
		else if (!*fname)
			fname = 0;
		if (fname)
			ADVari::debug_file = fopen(fname, "w");
		ADVari::zap_gcgen1 = -1;
		}
#endif
	if ((d = DErp::LastDerp) != 0) {
		for(i = 0; i < n; i++)
			V[i]->cv->aval = w[i];
#ifdef RAD_DEBUG
		if (ADVari::debug_file)
			do {
				fprintf(ADVari::debug_file, "%d\t%d\t%g + %g * %g",
					d->c->opno, d->b->opno, d->c->aval, *d->a, d->b->aval);
				d->c->aval += *d->a * d->b->aval;
				fprintf(ADVari::debug_file, " = %g\n", d->c->aval);
				fflush(ADVari::debug_file);
				} while((d = d->next));
		else
#endif
		do d->c->aval += *d->a * d->b->aval;
		while((d = d->next));
		}
#ifdef RAD_DEBUG
	if (ADVari::debug_file) {
		fclose(ADVari::debug_file);
		ADVari::debug_file = 0;
		}
#endif //RAD_DEBUG
#ifdef RAD_AUTO_AD_Const
	*ADVari::Last_ADvari = 0;
	ADVari::Last_ADvari = &ADVari::First_ADvari;
	if ((anext = ADVari::First_ADvari) && !(ADVari::adc.rad_need_reinit & 2)) {
		ADVari::adc.rad_need_reinit = 3;
		while((a = anext)) {
			anext = a->Next;
			if ((v = (IndepADvar<Double> *)a->padv)) {
#ifdef RAD_Const_WARN
				if ((i = a->opno) > 0)
					i = -i;
				j = a->gcgen;
				v->cv = cv = new ADVari(v, a->Val, a->aval);
				cv->opno = i;
				cv->gcgen = j;
#else
				v->cv = new ADVari(v, a->Val, a->aval);
#endif
				}
			}
		DErp::LastDerp = 0;
		}
#endif /*RAD_AUTO_AD_Const*/
#ifdef RAD_DEBUG
	ADVari::gcgen_cur++;
	ADVari::last_opno = 0;
#endif
		}

 template<typename Double>
IndepADvar<Double>::IndepADvar(Ttype d)
{

	ADVari *x = new ADVari(Hv_const, d);
	cv = x;
	}

 template<typename Double>
IndepADvar<Double>::IndepADvar(double i)
{

	ADVari *x = new ADVari(Hv_const, Double(i));
	cv = x;
	}

 template<typename Double>
IndepADvar<Double>::IndepADvar(int i)
{

	ADVari *x = new ADVari(Hv_const, Double(i));
	cv = x;
	}

 template<typename Double>
IndepADvar<Double>::IndepADvar(long i)
{

	ADVari *x = new ADVari(Hv_const, Double(i));
	cv = x;
	}

 template<typename Double>
ConstADvar<Double>::ConstADvar()
{
	ConstADVari *x = new ConstADVari(0.);
	this->cv = x;
	}

 template<typename Double> void
ConstADvar<Double>::ConstADvar_ctr(Double d)
{
	ConstADVari *x = new ConstADVari(d);
	this->cv = x;
	}

 template<typename Double>
ConstADvar<Double>::ConstADvar(const IndepADVar &x)
{
	ConstADVari *y = new ConstADVari(x.cv->Val);
	DErp *d = new DErp(&x.adc.One, y, x.cv);
	this->cv = y;
	}

 template<typename Double>
ConstADvar<Double>::ConstADvar(const ConstADvar &x)
{
	ConstADVari *y = new ConstADVari(x.cv->Val);
	DErp *d = new DErp(&x.cv->adc.One, y, (ADVari*)x.cv);
	this->cv = y;
	}

 template<typename Double>
ConstADvar<Double>::ConstADvar(const ADVari &x)
{
	ConstADVari *y = new ConstADVari(x.Val);
	DErp *d = new DErp(&x.adc.One, y, &x);
	this->cv = y;
	}

 template<typename Double>
 void
IndepADvar<Double>::AD_Const(const IndepADvar &v)
{
	typedef ConstADvari<Double> ConstADVari;

	ConstADVari *ncv = new ConstADVari(v.val());
#ifdef RAD_AUTO_AD_Const
	if (v.cv)
		v.cv->padv = 0;
#endif
	v.cv = ncv;
	}

 template<typename Double>
 void
ConstADvari<Double>::aval_reset()
{
	ConstADvari *x = ConstADvari::lastcad;
	while(x) {
		x->aval = 0;
		x = x->prevcad;
		}
	}

#ifdef RAD_AUTO_AD_Const

 template<typename Double>
ADvari<Double>::ADvari(const IndepADVar *x, Double d): Val(d), aval(0.)
{
	opclass = Hv_const;
	this->ADvari_padv(x);
	}

 template<typename Double>
ADvari<Double>::ADvari(const IndepADVar *x, Double d, Double g): Val(d), aval(g)
{
	opclass = Hv_const;
	this->ADvari_padv(x);
	}

 template<typename Double>
ADvar1<Double>::ADvar1(const IndepADVar *x, const IndepADVar &y):
	ADVari(Hv_copy, y.cv->Val), d((const Double*)&ADcontext<Double>::One, (ADVari*)this, y.cv)
{
	this->ADvari_padv(x);
	}

 template<typename Double>
ADvar1<Double>::ADvar1(const IndepADVar *x, const ADVari &y):
	ADVari(Hv_copy, y.Val), d((const Double*)&ADcontext<Double>::One, this, &y)
{
	this->ADvari_padv(x);
	}

#endif /* RAD_AUTO_AD_Const */

 template<typename Double>
 IndepADvar<Double>&
ADvar_operatoreq(IndepADvar<Double> *This, const ADvari<Double> &x)
{ This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x);
  return *(IndepADvar<Double>*) This; }

 template<typename Double>
 ADvar<Double>&
ADvar_operatoreq(ADvar<Double> *This, const ADvari<Double> &x)
{ This->cv = new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x); return *(ADvar<Double>*) This; }

 template<typename Double>
 IndepADvar<Double>&
IndepADvar<Double>::operator=(Double d)
{
#ifdef RAD_AUTO_AD_Const
	ADVari *ncv = new ADVari(this, d);
	if (this->cv)
		this->cv->padv = 0;
	this->cv = ncv;
#else
	this->cv = new ADVari(Hv_const, d);
#endif
	return *this;
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator=(Double d)
{
#ifdef RAD_AUTO_AD_Const
	ADVari *nv = new ADVari(this, d);
	if (this->cv)
		this->cv->padv = 0;
	this->cv = nv;
#else
	this->cv = ConstADVari::cadc.fpval_implies_const
		? new ConstADVari(d)
		: new ADVari(Hv_const, d);
#endif
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator-(const ADvari<Double> &T) {
	return *(new ADvar1<Double>(Hv_negate, -T.Val, &T.adc.negOne, &T));
	}

 template<typename Double>
 ADvari<Double>&
operator+(const ADvari<Double> &L, const ADvari<Double> &R) {
	return *(new ADvar2<Double>(Hv_plusLR, L.Val + R.Val, &L, &L.adc.One, &R, &L.adc.One));
	}

#ifdef RAD_AUTO_AD_Const
#define RAD_ACA ,this
#else
#define RAD_ACA /*nothing*/
#endif

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator+=(const ADVari &R) {
	ADVari *Lcv = this->cv;
	this->cv = new ADvar2<Double>(Hv_plusLR, Lcv->Val + R.Val, Lcv,
					&R.adc.One, &R, &R.adc.One RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator+(const ADvari<Double> &L, Double R) {
	return *(new ADvar1<Double>(Hv_copy, L.Val + R, &L.adc.One, &L));
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator+=(Double R) {
	ADVari *tcv = this->cv;
	this->cv = new ADVar1(Hv_copy, tcv->Val + R, &tcv->adc.One, tcv RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator+(Double L, const ADvari<Double> &R) {
	return *(new ADvar1<Double>(Hv_copy, L + R.Val, &R.adc.One, &R));
	}

 template<typename Double>
 ADvari<Double>&
operator-(const ADvari<Double> &L, const ADvari<Double> &R) {
	return *(new ADvar2<Double>(Hv_minusLR, L.Val - R.Val, &L, &L.adc.One, &R, &L.adc.negOne));
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator-=(const ADVari &R) {
	ADVari *Lcv = this->cv;
	this->cv = new ADvar2<Double>(Hv_minusLR,Lcv->Val - R.Val, Lcv,
					&R.adc.One, &R, &R.adc.negOne RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator-(const ADvari<Double> &L, Double R) {
	return *(new ADvar1<Double>(Hv_copy, L.Val - R, &L.adc.One, &L));
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator-=(Double R) {
	ADVari *tcv = this->cv;
	this->cv = new ADVar1(Hv_copy, tcv->Val - R, &tcv->adc.One, tcv RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator-(Double L, const ADvari<Double> &R) {
	return *(new ADvar1<Double>(Hv_negate, L - R.Val, &R.adc.negOne, &R));
	}

 template<typename Double>
 ADvari<Double>&
operator*(const ADvari<Double> &L, const ADvari<Double> &R) {
	return *(new ADvar2<Double>(Hv_timesLR, L.Val * R.Val, &L, &R.Val, &R, &L.Val));
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator*=(const ADVari &R) {
	ADVari *Lcv = this->cv;
	this->cv = new ADvar2<Double>(Hv_timesLR, Lcv->Val * R.Val, Lcv,
					&R.Val, &R, &Lcv->Val RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator*(const ADvari<Double> &L, Double R) {
	return *(new ADvar1s<Double>(L.Val * R, R, &L));
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator*=(Double R) {
	ADVari *Lcv = this->cv;
	this->cv = new ADvar1s<Double>(Lcv->Val * R, R, Lcv RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator*(Double L, const ADvari<Double> &R) {
	return *(new ADvar1s<Double>(L * R.Val, L, &R));
	}

 template<typename Double>
 ADvari<Double>&
operator/(const ADvari<Double> &L, const ADvari<Double> &R) {
	Double Lv = L.Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
	return *(new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, &L, &R));
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator/=(const ADVari &R) {
	ADVari *Lcv = this->cv;
	Double Lv = Lcv->Val, Rv = R.Val, pL = 1. / Rv, q = Lv/Rv, qpL = q*pL;
	this->cv = new ADvar2q<Double>(q, pL, -qpL, -pL*pL, 2.*pL*qpL, Lcv, &R RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
operator/(const ADvari<Double> &L, Double R) {
	return *(new ADvar1s<Double>(L.Val / R, 1./R, &L));
	}

 template<typename Double>
 ADvari<Double>&
operator/(Double L, const ADvari<Double> &R) {
	Double recip = 1. / R.Val;
	Double q = L * recip;
	Double d1 = -q*recip;
	return *(new ADvar1g<Double>(q, d1, -q*d1, &R));
	}

 template<typename Double>
 ADvar<Double>&
ADvar<Double>::operator/=(Double R) {
	ADVari *Lcv = this->cv;
	this->cv = new ADvar1s<Double>(Lcv->Val / R, 1./R, Lcv RAD_ACA);
	return *this;
	}

 template<typename Double>
 ADvari<Double>&
acos(const ADvari<Double> &v) {
	Double t = v.Val, t1 = 1. - t*t, d1 = -1./std::sqrt(t1);
	return *(new ADvar1g<Double>(std::acos(t), d1, t*d1/t1, &v));
	}

 template<typename Double>
 ADvari<Double>&
acosh(const ADvari<Double> &v) {
	Double d1, t, t1, t2;
	t = v.Val;
	t1 = std::sqrt(t2 = t*t - 1.);
	d1 = 1. / t1;
	return *(new ADvar1g<Double>(std::log(t + t1), d1, -t*d1/t2, &v));
	}

 template<typename Double>
 ADvari<Double>&
asin(const ADvari<Double> &v) {
	Double d1, t, t1;
	t = v.Val;
	d1 = 1. / std::sqrt(t1 = 1. - t*t);
	return *(new ADvar1g<Double>(std::asin(t), d1, t*d1/t1, &v));
	}

 template<typename Double>
 ADvari<Double>&
asinh(const ADvari<Double> &v) {
	Double d1, t, t1, t2, td;
	t = v.Val;
	t1 = std::sqrt(t2 = t*t + 1.);
	d1 = 1. / t1;
	td = 1.;
	if (t < 0.)
		td = -1.;
	return *(new ADvar1g<Double>(td*std::log(t*td + t1), d1, -(t/t2)*d1, &v));
	}

 template<typename Double>
 ADvari<Double>&
atan(const ADvari<Double> &v) {
	Double t = v.Val, d1 = 1./(1. + t*t);
	return *(new ADvar1g<Double>(std::atan(t), d1, -(t+t)*d1*d1, &v));
	}

 template<typename Double>
 ADvari<Double>&
atanh(const ADvari<Double> &v) {
	Double t = v.Val, d1 = 1./(1. - t*t);
	return *(new ADvar1g<Double>(0.5*std::log((1.+t)/(1.-t)), d1, (t+t)*d1*d1, &v));
	}

 template<typename Double>
 ADvari<Double>&
atan2(const ADvari<Double> &L, const ADvari<Double> &R) {
	Double R2, t, t2, x, x2, y, y2;
	x = L.Val;
	y = R.Val;
	x2 = x*x;
	y2 = y*y;
	t = 1./(x2 + y2);
	t2 = t*t;
	R2 = 2.*t2*x*y;
	return *(new ADvar2g<Double>(std::atan2(x,y), y*t, -x*t, -R2, t2*(x2 - y2), R2, &L, &R));
	}

 template<typename Double>
 ADvari<Double>&
atan2(Double x, const ADvari<Double> &R) {
	Double t, x2, y, y2;
	y = R.Val;
	x2 = x*x;
	y2 = y*y;
	t = 1./(x2 + y2);
	return *(new ADvar1g<Double>(std::atan2(x,y), -x*t, 2.*t*t*x*y, &R));
	}

 template<typename Double>
 ADvari<Double>&
atan2(const ADvari<Double> &L, Double y) {
	Double t, x, x2, y2;
	x = L.Val;
	x2 = x*x;
	y2 = y*y;
	t = 1./(x2 + y2);
	return *(new ADvar1g<Double>(std::atan2(x,y), y*t, -2.*t*t*x*y, &L));
	}

 template<typename Double>
 ADvari<Double>&
max(const ADvari<Double> &L, const ADvari<Double> &R) {
	const ADvari<Double> &x = L.Val >= R.Val ? L : R;
	return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
	}

 template<typename Double>
 ADvari<Double>&
max(Double L, const ADvari<Double> &R) {
	if (L >= R.Val)
		return *(new ADvari<Double>(Hv_const, L));
	return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
	}

 template<typename Double>
 ADvari<Double>&
max(const ADvari<Double> &L, Double R) {
	if (L.Val >= R)
		return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
	return *(new ADvari<Double>(Hv_const, R));
	}

 template<typename Double>
 ADvari<Double>&
min(const ADvari<Double> &L, const ADvari<Double> &R) {
	const ADvari<Double> &x = L.Val <= R.Val ? L : R;
	return *(new ADvar1<Double>(Hv_copy, x.Val, &x.adc.One, &x));
	}

 template<typename Double>
 ADvari<Double>&
min(Double L, const ADvari<Double> &R) {
	if (L <= R.Val)
		return *(new ADvari<Double>(Hv_const, L));
	return *(new ADvar1<Double>(Hv_copy, R.Val, &R.adc.One, &R));
	}

 template<typename Double>
 ADvari<Double>&
min(const ADvari<Double> &L, Double R) {
	if (L.Val <= R)
		return *(new ADvar1<Double>(Hv_copy, L.Val, &L.adc.One, &L));
	return *(new ADvari<Double>(Hv_const, R));
	}

 template<typename Double>
 ADvari<Double>&
cos(const ADvari<Double> &v) {
	Double t = std::cos(v.Val);
	return *(new ADvar1g<Double>(t, -std::sin(v.Val), -t, &v));
	}

 template<typename Double>
 ADvari<Double>&
cosh(const ADvari<Double> &v) {
	Double t = std::cosh(v.Val);
	return *(new ADvar1g<Double>(t, std::sinh(v.Val), t, &v));
	}

 template<typename Double>
 ADvari<Double>&
exp(const ADvari<Double> &v) {
	Double t = std::exp(v.Val);
	return *(new ADvar1g<Double>(t, t, t, &v));
	}

 template<typename Double>
 ADvari<Double>&
log(const ADvari<Double> &v) {
	Double x = v.Val, d1 = 1. / x;
	return *(new ADvar1g<Double>(std::log(x), d1, -d1*d1, &v));
	}

 template<typename Double>
 ADvari<Double>&
log10(const ADvari<Double> &v) {
	static double num = 1. / std::log(10.);
	Double d1, t, x;
	x = v.Val;
	t = 1. / x;
	d1 = num * t;
	return *(new ADvar1g<Double>(std::log10(x), d1, -d1*t, &v));
	}

 template<typename Double>
 ADvari<Double>&
pow(const ADvari<Double> &L, const ADvari<Double> &R) {
	Double dx, dy, t, x, xlog, xym1, y;
	x = L.Val;
	y = R.Val;
	t = std::pow(x,y);
	xym1 = t / x;
	xlog = std::log(x);
	dx = y*xym1;
	dy = t * xlog;
	return *(new ADvar2g<Double>(t, dx, dy, (y-1.)*dx/x, xym1*(1. + y*xlog), dy*xlog, &L, &R));
	}

 template<typename Double>
 ADvari<Double>&
pow(Double x, const ADvari<Double> &R) {
	Double dy, t, xlog, y;
	y = R.Val;
	t = std::pow(x,y);
	xlog = std::log(x);
	dy = t * xlog;
	return *(new ADvar1g<Double>(t, dy, dy*xlog, &R));
	}

 template<typename Double>
 ADvari<Double>&
pow(const ADvari<Double> &L, Double y) {
	Double dx, t, x;
	x = L.Val;
	t = std::pow(x,y);
	dx = y*t/x;
	return *(new ADvar1g<Double>(t, dx, (y-1.)*dx/x, &L));
	}

 template<typename Double>
 ADvari<Double>&
sin(const ADvari<Double> &v) {
	Double t = std::sin(v.Val);
	return *(new ADvar1g<Double>(t, std::cos(v.Val), -t, &v));
	}

 template<typename Double>
 ADvari<Double>&
sinh(const ADvari<Double> &v) {
	Double t = std::sinh(v.Val);
	return *(new ADvar1g<Double>(t, std::cosh(v.Val), t, &v));
	}

 template<typename Double>
 ADvari<Double>&
sqrt(const ADvari<Double> &v) {
	Double t = std::sqrt(v.Val);
	Double d1 = 0.5 / t;
	return *(new ADvar1g<Double>(t, d1, -0.5*d1/v.Val, &v));
	}

 template<typename Double>
 ADvari<Double>&
tan(const ADvari<Double> &v) {
	Double d1, rv, t;
	rv = std::tan(v.Val);
	t = 1. / std::cos(v.Val);
	d1 = t*t;
	return *(new ADvar1g<Double>(rv, d1, (rv+rv)*d1, &v));
	}

 template<typename Double>
 ADvari<Double>&
tanh(const ADvari<Double> &v) {
	Double d1, rv, t;
	rv = std::tanh(v.Val);
	t = 1. / std::cosh(v.Val);
	d1 = t*t;
	return *(new ADvar1g<Double>(rv, d1, -(rv+rv)*d1, &v));
	}

 template<typename Double>
 ADvari<Double>&
abs(const ADvari<Double> &v) {
	Double t, p;
	p = 1.;
	if ((t = v.Val) < 0) {
		t = -t;
		p = -p;
		}
	return *(new ADvar1g<Double>(t, p, 0., &v));
	}

 template<typename Double>
 ADvari<Double>&
fabs(const ADvari<Double> &v) {	// Synonym for "abs"
				// "fabs" is not the best choice of name,
				// but this name is used at Sandia.
	Double t, p;
	p = 1.;
	if ((t = v.Val) < 0) {
		t = -t;
		p = -p;
		}
	return *(new ADvar1g<Double>(t, p, 0., &v));
	}

 template<typename Double>
 ADvari<Double>&
ADf1(Double f, Double g, Double h, const ADvari<Double> &x) {
	return *(new ADvar1g<Double>(f, g, h, &x));
	}

 template<typename Double>
 inline ADvari<Double>&
ADf1(Double f, Double g, Double h, const IndepADvar<Double> &x) {
	return *(new ADvar1g<Double>(f, g, h, x.cv));
	}

 template<typename Double>
 ADvari<Double>&
ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
		const ADvari<Double> &x, const ADvari<Double> &y) {
	return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, &y));
	}

 template<typename Double>
 ADvari<Double>&
ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
		const ADvari<Double> &x, const IndepADvar<Double> &y) {
	return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, &x, y.cv));
	}

 template<typename Double>
 ADvari<Double>&
ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
		const IndepADvar<Double> &x, const ADvari<Double> &y) {
	return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, &y));
	}

 template<typename Double>
 ADvari<Double>&
ADf2(Double f, Double gx, Double gy, Double hxx, Double hxy, Double hyy,
		const IndepADvar<Double> &x, const IndepADvar<Double> &y) {
	return *(new ADvar2g<Double>(f, gx, gy, hxx, hxy, hyy, x.cv, y.cv));
	}

 template<typename Double>
 ADvari<Double>&
ADfn(Double f, int n, const IndepADvar<Double> *x, const Double *g, const Double *h) {
	return *(new ADvarn<Double>(f, n, x, g, h));
	}

 template<typename Double>
 inline ADvari<Double>&
ADfn(Double f, int n, const ADvar<Double> *x, const Double *g, const Double *h) {
	return ADfn<Double>(f, n, (IndepADvar<Double>*)x, g, h);
	}

 template<typename Double>
 void
ADcontext<Double>::Hvprod(int n, ADvar<Double> **x, Double *v, Double *hv)
{
	ADVari *a, *aL, *aR, **ap, **ape;
	ADVari_block *b, *b0;
	DErp *d;
	Double aO, adO, *g, *h, *h0, t, tL, tR;
	int i, j, k, m;
	for(i = 0; i < n; i++) {
		a = x[i]->cv;
		a->dO = v[i];
		a->aO = a->adO = 0.;
		}
	ADVari::adc.Aibusy->limit = ADVari::adc.Ainext;
	a = 0;
	for(b0 = 0, b = ADVari::adc.AiFirst; b; b0 = b, b = b->next) {
		ap = b->pADvari;
		ape = b->limit;
		while(ap < ape) {
			a = *ap++;
			a->aO = a->adO = 0.;
			switch(a->opclass) {
			 case Hv_copy:
				a->dO = ((ADVar1*)a)->d.c->dO;
				break;
			 case Hv_binary:
				a->dO =   ((ADVar2g*)a)->pL * ((ADVar2g*)a)->dL.c->dO
					+ ((ADVar2g*)a)->pR * ((ADVar2g*)a)->dR.c->dO;
				break;
			 case Hv_unary:
				a->dO = ((ADVar1g*)a)->pL * ((ADVar1g*)a)->d.c->dO;
				break;
			 case Hv_negate:
				a->dO = -((ADVar1*)a)->d.c->dO;
				break;
			 case Hv_plusLR:
				a->dO = ((ADVar2*)a)->dL.c->dO + ((ADVar2*)a)->dR.c->dO;
				break;
			 case Hv_minusLR:
				a->dO = ((ADVar2*)a)->dL.c->dO - ((ADVar2*)a)->dR.c->dO;
				break;
			 case Hv_timesL:
				a->dO = ((ADVar1s*)a)->pL * ((ADVar1s*)a)->d.c->dO;
				break;
			 case Hv_timesLR:
				a->dO =   ((ADVar2*)a)->dR.c->Val * ((ADVar2*)a)->dL.c->dO
					+ ((ADVar2*)a)->dL.c->Val * ((ADVar2*)a)->dR.c->dO;
				break;
			 case Hv_quotLR:
				a->dO =   ((ADVar2q*)a)->pL * ((ADVar2q*)a)->dL.c->dO
					+ ((ADVar2q*)a)->pR * ((ADVar2q*)a)->dR.c->dO;
				break;
			 case Hv_nary:
				d = ((ADVarn*)a)->D;
				m = ((ADVarn*)a)->n;
				g = ((ADVarn*)a)->G;
				t = 0.;
				for(i = 0; i < m; i++)
					t += g[i] * d[i].c->dO;
				a->dO = t;
                                break;
                        case Hv_const:
                          ;
			 }
			}
		}
	if (a)
		a->adO = 1.;
	for(b = b0; b; b = b->prev) {
		ape = b->pADvari;
		ap = b->limit;
		while(ap > ape) {
			a = *--ap;
			aO = a->aO;
			adO = a->adO;
			switch(a->opclass) {
			 case Hv_copy:
				aL = ((ADVar1*)a)->d.c;
				aL->aO += aO;
				aL->adO += adO;
				break;
			 case Hv_binary:
				aL = ((ADVar2g*)a)->dL.c;
				aR = ((ADVar2g*)a)->dR.c;
				tL = adO*aL->dO;
				tR = adO*aR->dO;
				aL->aO += aO*((ADVar2g*)a)->pL
					+ tL*((ADVar2g*)a)->pL2
					+ tR*((ADVar2g*)a)->pLR;
				aR->aO += aO*((ADVar2g*)a)->pR
					+ tL*((ADVar2g*)a)->pLR
					+ tR*((ADVar2g*)a)->pR2;
				aL->adO += adO * ((ADVar2g*)a)->pL;
				aR->adO += adO * ((ADVar2g*)a)->pR;
				break;
			 case Hv_unary:
				aL = ((ADVar1g*)a)->d.c;
				aL->aO += aO * ((ADVar1g*)a)->pL
					+ adO * aL->dO * ((ADVar1g*)a)->pL2;
				aL->adO += adO * ((ADVar1g*)a)->pL;
				break;
			 case Hv_negate:
				aL = ((ADVar1*)a)->d.c;
				aL->aO -= aO;
				aL->adO -= adO;
				break;
			 case Hv_plusLR:
				aL = ((ADVar2*)a)->dL.c;
				aR = ((ADVar2*)a)->dR.c;
				aL->aO += aO;
				aL->adO += adO;
				aR->aO += aO;
				aR->adO += adO;
				break;
			 case Hv_minusLR:
				aL = ((ADVar2*)a)->dL.c;
				aR = ((ADVar2*)a)->dR.c;
				aL->aO += aO;
				aL->adO += adO;
				aR->aO -= aO;
				aR->adO -= adO;
				break;
			 case Hv_timesL:
				aL = ((ADVar1s*)a)->d.c;
				aL->aO += aO * (tL = ((ADVar1s*)a)->pL);
				aL->adO += adO * tL;
				break;
			 case Hv_timesLR:
				aL = ((ADVar2*)a)->dL.c;
				aR = ((ADVar2*)a)->dR.c;
				aL->aO += aO * (tL = aR->Val) + adO*aR->dO;
				aR->aO += aO * (tR = aL->Val) + adO*aL->dO;
				aL->adO += adO * tL;
				aR->adO += adO * tR;
				break;
			 case Hv_quotLR:
				aL = ((ADVar2q*)a)->dL.c;
				aR = ((ADVar2q*)a)->dR.c;
				tL = adO*aL->dO;
				tR = adO*aR->dO;
				aL->aO += aO*((ADVar2q*)a)->pL
					+ tR*((ADVar2q*)a)->pLR;
				aR->aO += aO*((ADVar2q*)a)->pR
					+ tL*((ADVar2q*)a)->pLR
					+ tR*((ADVar2q*)a)->pR2;
				aL->adO += adO * ((ADVar2q*)a)->pL;
				aR->adO += adO * ((ADVar2q*)a)->pR;
				break;
			 case Hv_nary:
				d  = ((ADVarn*)a)->D;
				m  = ((ADVarn*)a)->n;
				g  = ((ADVarn*)a)->G;
				h0 = ((ADVarn*)a)->H;
				for(i = 0; i < m; i++) {
					aL = d[i].c;
					aL->adO += adO * (t = g[i]);
					aL->aO += t*aO;
					t = adO * aL->dO;
					for(h = h0, j = 0; j <= i; j++)
						d[j].c->aO += t * *h++;
					h0 = h--;
					for(k = j; j < m; j++)
						d[j].c->aO += t * *(h += k++);
					}
                        case Hv_const:
                          ;
			 }
			}
		}
	for(i = 0; i < n; i++) {
		a = x[i]->cv;
		a->dO = 0.;
		hv[i] = a->aO;
		}
	}

 template<typename Double>
 inline Double
val(const ADvari<Double> &x) {
	return x.Val;
	}

#undef RAD_ACA
#define A (ADvari<Double>*)
#ifdef RAD_Const_WARN
#define C(x) (((x)->opno < 0) ? RAD_Const_Warn(x) : 0, *A x)
#else
#define C(x) *A x
#endif
#define T template<typename Double> inline
#define F ADvari<Double>&
#define Ai const ADvari<Double>&
#define AI const IndepADvar<Double>&
#define D Double
#define T2(r,f) \
 T r f(Ai L, AI R) { return f(L, C(R.cv)); }\
 T r f(AI L, Ai R) { return f(C(L.cv), R); }\
 T r f(AI L, AI R) { return f(C(L.cv), C(R.cv)); }\
 T r f(AI L, D R) { return f(C(L.cv), R); }\
 T r f(Ai L, Dtype R) { return f(L, (D)R); }\
 T r f(AI L, Dtype R) { return f(C(L.cv), (D)R); }\
 T r f(Ai L, long R) { return f(L, (D)R); }\
 T r f(AI L, long R) { return f(C(L.cv), (D)R); }\
 T r f(Ai L, int R) { return f(L, (D)R); }\
 T r f(AI L, int R) { return f(C(L.cv), (D)R); }\
 T r f(D L, AI R) { return f(L, C(R.cv)); }\
 T r f(Dtype L, Ai R) { return f((D)L, R); }\
 T r f(Dtype L, AI R) { return f((D)L, C(R.cv)); }\
 T r f(long L, Ai R) { return f((D)L, R); }\
 T r f(long L, AI R) { return f((D)L, C(R.cv)); }\
 T r f(int L, Ai R) { return f((D)L, R); }\
 T r f(int L, AI R) { return f((D)L, C(R.cv)); }

T2(F, operator+)
T2(F, operator-)
T2(F, operator*)
T2(F, operator/)
T2(F, atan2)
T2(F, pow)
T2(F, max)
T2(F, min)
T2(int, operator<)
T2(int, operator<=)
T2(int, operator==)
T2(int, operator!=)
T2(int, operator>=)
T2(int, operator>)

#undef T2
#undef D

#define T1(f)\
 T F f(AI x) { return f(C(x.cv)); }

T1(operator+)
T1(operator-)
T1(abs)
T1(acos)
T1(acosh)
T1(asin)
T1(asinh)
T1(atan)
T1(atanh)
T1(cos)
T1(cosh)
T1(exp)
T1(log)
T1(log10)
T1(sin)
T1(sinh)
T1(sqrt)
T1(tan)
T1(tanh)
T1(fabs)

T F copy(AI x)
{ return *(new ADvar1<Double>(Hv_copy, x.cv->Val, &ADcontext<Double>::One, (ADvari<Double>*)x.cv)); }

T F copy(Ai x)
{ return *(new ADvar1<Double>(Hv_copy, x.Val, &ADcontext<Double>::One, (ADvari<Double>*)&x)); }

#undef T1
#undef AI
#undef Ai
#undef F
#undef T
#undef A
#undef C
#undef Ttype
#undef Dtype

} /* namespace Rad2 */
} /* namespace Sacado */

#endif /* SACADO_TRAD2_H */
