// Copyright (C) 1996-2007 Ole Stauning & Claus Bendtsen (fadbad@uning.dk) // All rights reserved. // This code is provided "as is", without any warranty of any kind, // either expressed or implied, including but not limited to, any implied // warranty of merchantibility or fitness for any purpose. In no event // will any party who distributed the code be liable for damages or for // any claim(s) by any other party, including but not limited to, any // lost profits, lost monies, lost data or data rendered inaccurate, // losses sustained by third parties, or any other special, incidental or // consequential damages arising out of the use or inability to use the // program, even if the possibility of such damages has been advised // against. The entire risk as to the quality, the performance, and the // fitness of the program for any particular purpose lies with the party // using the code. // This code, and any derivative of this code, may not be used in a // commercial package without the prior explicit written permission of // the authors. Verbatim copies of this code may be made and distributed // in any medium, provided that this copyright notice is not removed or // altered in any way. No fees may be charged for distribution of the // codes, other than a fee to cover the cost of the media and a // reasonable handling fee. // *************************************************************** // ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE // COPYRIGHT NOTICE // *************************************************************** #ifndef _FADIFF_H #define _FADIFF_H #include "fadbad.h" namespace fadbad { template class FTypeName // STACK-BASED { T m_val; T m_diff[N]; bool m_depend; public: typedef T UnderlyingType; FTypeName():m_depend(false){} FTypeName(const FTypeName& val):m_val(val.m_val),m_depend(val.m_depend) { if (m_depend) for(unsigned int i=0;i /*explicit*/ FTypeName(const U& val):m_val(val),m_depend(false) { } template FTypeName& operator=(const U& val) { m_val=val; m_depend=false; return *this; } FTypeName& operator=(const FTypeName& val) { if (this==&val) return *this; m_val=val.m_val; m_depend=val.m_depend; if (m_depend) for(unsigned int i=0;i::myZero(); return zero; } T& d(const unsigned int i) { USER_ASSERT(i::myZero(); return zero; } T& diff(unsigned int idx) { USER_ASSERT(idx::myZero(); m_diff[i++]=Op::myOne(); for( ;i::myZero(); m_depend=true; return m_diff[idx]; } bool depend() const { return m_depend; } void setDepend(const FTypeName&) { m_depend=true; } void setDepend(const FTypeName&, const FTypeName&) { m_depend=true; } FTypeName& operator+=(const FTypeName& val); FTypeName& operator-=(const FTypeName& val); FTypeName& operator*=(const FTypeName& val); FTypeName& operator/=(const FTypeName& val); template FTypeName& operator+=(const V& val); template FTypeName& operator-=(const V& val); template FTypeName& operator*=(const V& val); template FTypeName& operator/=(const V& val); }; template class FTypeName // HEAP-BASED { T m_val; unsigned int m_size; T* m_diff; public: typedef T UnderlyingType; FTypeName():m_val(),m_size(0),m_diff(0){} FTypeName(const FTypeName& val):m_val(val.m_val),m_size(val.m_size),m_diff(m_size==0?0:new T[m_size]) { for(unsigned int i=0;i /*explicit*/ FTypeName(const U& val):m_val(val),m_size(0),m_diff(0){} ~FTypeName(){ delete[] m_diff; } template FTypeName& operator=(const U& val) { m_val=val; m_size=0; delete[] m_diff; m_diff=0; return *this; } FTypeName& operator=(const FTypeName& val) { if (this==&val) return *this; m_val=val.m_val; if (val.m_size>0) { if (m_size==0) { m_size = val.m_size; m_diff = new T[m_size]; } USER_ASSERT(m_size==val.m_size,"derivative vectors not of same size"); for(unsigned int i=0;i0) { for(unsigned int i=0;i::myZero(); return zero; } T& d(const unsigned int i) { if (i::myZero(); return zero; } T& diff(unsigned int idx, unsigned int N) { USER_ASSERT(idx::myZero(); m_diff[i++]=Op::myOne(); for( ;i::myZero(); return m_diff[idx]; } bool depend() const { return m_size!=0; } void setDepend(const FTypeName& val) { INTERNAL_ASSERT(val.m_size>0,"input is not a dependent variable") if (m_size==0) { m_size = val.m_size; m_diff = new T[m_size]; } else { USER_ASSERT(m_size==val.m_size,"derivative vectors not of same size "<& val1, const FTypeName& val2) { USER_ASSERT(val1.m_size==val2.m_size,"derivative vectors not of same size "<0,"lhs-input is not a dependent variable") INTERNAL_ASSERT(val2.m_size>0,"rhs-input is not a dependent variable") if (m_size==0) { m_size=val1.m_size; m_diff = new T[m_size]; } else { USER_ASSERT(m_size==val1.m_size,"derivative vectors not of same size "<& operator+=(const FTypeName& val); FTypeName& operator-=(const FTypeName& val); FTypeName& operator*=(const FTypeName& val); FTypeName& operator/=(const FTypeName& val); template FTypeName& operator+=(const V& val); template FTypeName& operator-=(const V& val); template FTypeName& operator*=(const V& val); template FTypeName& operator/=(const V& val); }; template bool operator==(const FTypeName& val1, const FTypeName& val2) { return Op::myEq(val1.val(),val2.val()); } template bool operator!=(const FTypeName& val1, const FTypeName& val2) { return Op::myNe(val1.val(),val2.val()); } template bool operator<(const FTypeName& val1, const FTypeName& val2) { return Op::myLt(val1.val(),val2.val()); } template bool operator<=(const FTypeName& val1, const FTypeName& val2) { return Op::myLe(val1.val(),val2.val()); } template bool operator>(const FTypeName& val1, const FTypeName& val2) { return Op::myGt(val1.val(),val2.val()); } template bool operator>=(const FTypeName& val1, const FTypeName& val2) { return Op::myGe(val1.val(),val2.val()); } template bool operator==(const FTypeName& val1, const U& val2) { return Op::myEq(val1.val(),val2); } template bool operator==(const U& val1, const FTypeName& val2) { return Op::myEq(val1,val2.val()); } template bool operator!=(const FTypeName& val1, const U& val2) { return Op::myNe(val1.val(),val2); } template bool operator!=(const U& val1, const FTypeName& val2) { return Op::myNe(val1,val2.val()); } template bool operator<(const FTypeName& val1, const U& val2) { return Op::myLt(val1.val(),val2); } template bool operator<(const U& val1, const FTypeName& val2) { return Op::myLt(val1,val2.val()); } template bool operator<=(const FTypeName& val1, const U& val2) { return Op::myLe(val1.val(),val2); } template bool operator<=(const U& val1, const FTypeName& val2) { return Op::myLe(val1,val2.val()); } template bool operator>(const FTypeName& val1, const U& val2) { return Op::myGt(val1.val(),val2); } template bool operator>(const U& val1, const FTypeName& val2) { return Op::myGt(val1,val2.val()); } template bool operator>=(const FTypeName& val1, const U& val2) { return Op::myGe(val1.val(),val2); } template bool operator>=(const U& val1, const FTypeName& val2) { return Op::myGe(val1,val2.val()); } template INLINE2 FTypeName add1(const U& a, const FTypeName& b) { FTypeName c(a+b.val()); if (!b.depend()) return c; c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName add1(const U& a, const FTypeName& b) { FTypeName c(a+b.val()); if (!b.depend()) return c; c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName add2(const FTypeName& a, const U& b) { FTypeName c(a.val()+b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName add2(const FTypeName& a, const U& b) { FTypeName c(a.val()+b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName operator+ (const U& a, const FTypeName& b) { return add1(a,b); } template INLINE2 FTypeName operator+ (const FTypeName& a, const U& b) { return add2(a,b); } template INLINE2 FTypeName add3(const FTypeName& a, const FTypeName& b) { FTypeName c(a.val()+b.val()); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName add3(const FTypeName& a, const FTypeName& b) { FTypeName c(a.val()+b.val()); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName operator+ (const FTypeName& a, const FTypeName& b) { switch ((a.depend()?1:0)|(b.depend()?2:0)) { case 0: return FTypeName(a.val()+b.val()); case 1: return add2(a,b.val()); case 2: return add1(a.val(),b); } return add3(a,b); } template INLINE2 FTypeName sub1(const U& a, const FTypeName& b) { FTypeName c(a-b.val()); if (!b.depend()) return c; c.setDepend(b); for(unsigned int i=0;i::myNeg(b[i]); return c; } template INLINE2 FTypeName sub1(const U& a, const FTypeName& b) { FTypeName c(a-b.val()); if (!b.depend()) return c; c.setDepend(b); for(unsigned int i=0;i::myNeg(b[i]); return c; } template INLINE2 FTypeName sub2(const FTypeName& a, const U& b) { FTypeName c(a.val()-b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName sub2(const FTypeName& a, const U& b) { FTypeName c(a.val()-b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName operator- (const U& a, const FTypeName& b) { return sub1(a,b); } template INLINE2 FTypeName operator- (const FTypeName& a, const U& b) { return sub2(a,b); } template INLINE2 FTypeName sub3(const FTypeName& a, const FTypeName& b) { FTypeName c(a.val()-b.val()); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName sub3(const FTypeName& a, const FTypeName& b) { FTypeName c(a.val()-b.val()); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName operator- (const FTypeName& a, const FTypeName& b) { switch ((a.depend()?1:0)|(b.depend()?2:0)) { case 0: return FTypeName(a.val()-b.val()); case 1: return sub2(a,b.val()); case 2: return sub1(a.val(),b); } return sub3(a,b); } template INLINE2 FTypeName mul1(const U& a, const FTypeName& b) { FTypeName c(a*b.val()); if (!b.depend()) return c; c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName mul1(const U& a, const FTypeName& b) { FTypeName c(a*b.val()); if (!b.depend()) return c; c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName mul2(const FTypeName& a, const U& b) { FTypeName c(a.val()*b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName mul2(const FTypeName& a, const U& b) { FTypeName c(a.val()*b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName operator* (const U& a, const FTypeName& b) { return mul1(a,b); } template INLINE2 FTypeName operator* (const FTypeName& a, const U& b) { return mul2(a,b); } template INLINE2 FTypeName mul3 (const FTypeName& a, const FTypeName& b) { const T& aval(a.val()); const T& bval(b.val()); FTypeName c(aval*bval); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName mul3 (const FTypeName& a, const FTypeName& b) { const T& aval(a.val()); const T& bval(b.val()); FTypeName c(aval*bval); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName operator* (const FTypeName& a, const FTypeName& b) { switch ((a.depend()?1:0)|(b.depend()?2:0)) { case 0: return FTypeName(a.val()*b.val()); case 1: return mul2(a,b.val()); case 2: return mul1(a.val(),b); } return mul3(a,b); } template INLINE2 FTypeName div1(const U& a, const FTypeName& b) { FTypeName c(a/b.val()); if (!b.depend()) return c; T tmp(Op::myNeg(c.val()/b.val())); c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName div1(const U& a, const FTypeName& b) { FTypeName c(a/b.val()); if (!b.depend()) return c; T tmp(Op::myNeg(c.val()/b.val())); c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName div2(const FTypeName& a, const U& b) { FTypeName c(a.val()/b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName div2(const FTypeName& a, const U& b) { FTypeName c(a.val()/b); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName operator/ (const U& a, const FTypeName& b) { return div1(a,b); } template INLINE2 FTypeName operator/ (const FTypeName& a, const U& b) { return div2(a,b); } template INLINE2 FTypeName div3(const FTypeName& a, const FTypeName& b) { const T& bval(b.val()); FTypeName c(a.val()/bval); c.setDepend(a,b); const T& cval(c.val()); for(unsigned int i=0;i INLINE2 FTypeName div3(const FTypeName& a, const FTypeName& b) { const T& bval(b.val()); FTypeName c(a.val()/bval); c.setDepend(a,b); const T& cval(c.val()); for(unsigned int i=0;i INLINE2 FTypeName operator/ (const FTypeName& a, const FTypeName& b) { switch ((a.depend()?1:0)|(b.depend()?2:0)) { case 0: return FTypeName(a.val()/b.val()); case 1: return div2(a,b.val()); case 2: return div1(a.val(),b); } return div3(a,b); } template FTypeName& FTypeName::operator+=(const FTypeName& val) { Op::myCadd(m_val,val.m_val); if (!val.depend()) return *this; if (this->depend()) { for(unsigned int i=0;i::myCadd(m_diff[i],val[i]); } else { this->setDepend(val); for(unsigned int i=0;i FTypeName& FTypeName::operator-=(const FTypeName& val) { Op::myCsub(m_val,val.m_val); if (!val.depend()) return *this; if (this->depend()) { for(unsigned int i=0;i::myCsub(m_diff[i],val[i]); } else { this->setDepend(val); for(unsigned int i=0;i::myNeg(val[i]); } return *this; } template FTypeName& FTypeName::operator*=(const FTypeName& val) { if (this->depend() && val.depend()) { for(unsigned int i=0;idepend()) { for(unsigned int i=0;i::myCmul(m_diff[i],val.m_val); } else // (val.depend()) { this->setDepend(val); for(unsigned int i=0;i::myCmul(m_val,val.m_val); return *this; } template FTypeName& FTypeName::operator/=(const FTypeName& val) { Op::myCdiv(m_val,val.m_val); if (this->depend() && val.depend()) { for(unsigned int i=0;idepend()) { for(unsigned int i=0;i::myCdiv(m_diff[i],val.m_val); } else // (val.depend()) { this->setDepend(val); for(unsigned int i=0;i::myNeg(m_val*val.m_diff[i]/val.m_val); } return *this; } template template FTypeName& FTypeName::operator+=(const V& val) { Op::myCadd(m_val,val); return *this; } template template FTypeName& FTypeName::operator-=(const V& val) { Op::myCsub(m_val,val); return *this; } template template FTypeName& FTypeName::operator*=(const V& val) { Op::myCmul(m_val,val); if (!this->depend()) return *this; for(unsigned int i=0;i::myCmul(m_diff[i],val); return *this; } template template FTypeName& FTypeName::operator/=(const V& val) { Op::myCdiv(m_val,val); if (!this->depend()) return *this; for(unsigned int i=0;i::myCdiv(m_diff[i],val); return *this; } template FTypeName& FTypeName::operator+=(const FTypeName& val) { Op::myCadd(m_val,val.m_val); if (!val.depend()) return *this; if (this->depend()) { for(unsigned int i=0;isize();++i) Op::myCadd(m_diff[i],val[i]); } else { this->setDepend(val); for(unsigned int i=0;isize();++i) m_diff[i]=val[i]; } return *this; } template FTypeName& FTypeName::operator-=(const FTypeName& val) { Op::myCsub(m_val,val.m_val); if (!val.depend()) return *this; if (this->depend()) { for(unsigned int i=0;isize();++i) Op::myCsub(m_diff[i],val[i]); } else { this->setDepend(val); for(unsigned int i=0;isize();++i) m_diff[i]=Op::myNeg(val[i]); } return *this; } template FTypeName& FTypeName::operator*=(const FTypeName& val) { if (this->depend() && val.depend()) { for(unsigned int i=0;isize();++i) m_diff[i]=m_diff[i]*val.m_val+val.m_diff[i]*m_val; } else if (this->depend()) { for(unsigned int i=0;isize();++i) Op::myCmul(m_diff[i],val.m_val); } else // (val.depend()) { this->setDepend(val); for(unsigned int i=0;isize();++i) m_diff[i]=val.m_diff[i]*m_val; } Op::myCmul(m_val,val.m_val); return *this; } template FTypeName& FTypeName::operator/=(const FTypeName& val) { Op::myCdiv(m_val,val.m_val); if (this->depend() && val.depend()) { for(unsigned int i=0;isize();++i) m_diff[i]=(m_diff[i]-m_val*val.m_diff[i])/val.m_val; } else if (this->depend()) { for(unsigned int i=0;isize();++i) Op::myCdiv(m_diff[i],val.m_val); } else // (val.depend()) { this->setDepend(val); for(unsigned int i=0;isize();++i) m_diff[i]=Op::myNeg(m_val*val.m_diff[i]/val.m_val); } return *this; } template template FTypeName& FTypeName::operator+=(const V& val) { Op::myCadd(m_val,val); return *this; } template template FTypeName& FTypeName::operator-=(const V& val) { Op::myCsub(m_val,val); return *this; } template template FTypeName& FTypeName::operator*=(const V& val) { Op::myCmul(m_val,val); if (!this->depend()) return *this; for(unsigned int i=0;isize();++i) Op::myCmul(m_diff[i],val); return *this; } template template FTypeName& FTypeName::operator/=(const V& val) { Op::myCdiv(m_val,val); if (!this->depend()) return *this; for(unsigned int i=0;isize();++i) Op::myCdiv(m_diff[i],val); return *this; } template INLINE2 FTypeName pow1(const U& a, const FTypeName& b) { FTypeName c(Op::myPow(a,b.val())); if (!b.depend()) return c; T tmp(c.val()*Op::myLog(a)); c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName pow1(const U& a, const FTypeName& b) { FTypeName c(Op::myPow(a,b.val())); if (!b.depend()) return c; T tmp(c.val()*Op::myLog(a)); c.setDepend(b); for(unsigned int i=0;i INLINE2 FTypeName pow2(const FTypeName& a, const U& b) { FTypeName c(Op::myPow(a.val(),b)); if (!a.depend()) return c; T tmp(b*Op::myPow(a.val(),b-Op::myOne())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName pow2(const FTypeName& a, const U& b) { FTypeName c(Op::myPow(a.val(),b)); if (!a.depend()) return c; T tmp(b*Op::myPow(a.val(),b-Op::myOne())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName pow (const U& a, const FTypeName& b) { return pow1(a,b); } template INLINE2 FTypeName pow (const FTypeName& a,const U& b) { return pow2(a,b); } template INLINE2 FTypeName pow3(const FTypeName& a, const FTypeName& b) { FTypeName c(Op::myPow(a.val(),b.val())); T tmp(b.val()*Op::myPow(a.val(),b.val()-Op::myOne())),tmp1(c.val()*Op::myLog(a.val())); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName pow3(const FTypeName& a, const FTypeName& b) { FTypeName c(Op::myPow(a.val(),b.val())); T tmp(b.val()*Op::myPow(a.val(),b.val()-Op::myOne())),tmp1(c.val()*Op::myLog(a.val())); c.setDepend(a,b); for(unsigned int i=0;i INLINE2 FTypeName pow (const FTypeName& a, const FTypeName& b) { switch ((a.depend()?1:0)|(b.depend()?2:0)) { case 0: return FTypeName(Op::myPow(a.val(),b.val())); case 1: return pow2(a,b.val()); case 2: return pow1(a.val(),b); } return pow3(a,b); } /* Unary operators */ template INLINE2 FTypeName operator+ (const FTypeName& a) { FTypeName c(a.val()); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName operator+ (const FTypeName& a) { FTypeName c(a.val()); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName operator- (const FTypeName& a) { FTypeName c(Op::myNeg(a.val())); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i::myNeg(a[i]); return c; } template INLINE2 FTypeName operator- (const FTypeName& a) { FTypeName c(Op::myNeg(a.val())); if (!a.depend()) return c; c.setDepend(a); for(unsigned int i=0;i::myNeg(a[i]); return c; } template INLINE2 FTypeName sqr (const FTypeName& a) { FTypeName c(Op::mySqr(a.val())); if (!a.depend()) return c; T tmp(Op::myTwo()*a.val()); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName sqr (const FTypeName& a) { FTypeName c(Op::mySqr(a.val())); if (!a.depend()) return c; T tmp(Op::myTwo()*a.val()); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName exp (const FTypeName& a) { FTypeName c(Op::myExp(a.val())); if (!a.depend()) return c; c.setDepend(a); const T& cval(c.val()); for(unsigned int i=0;i INLINE2 FTypeName exp (const FTypeName& a) { FTypeName c(Op::myExp(a.val())); if (!a.depend()) return c; c.setDepend(a); const T& cval(c.val()); for(unsigned int i=0;i INLINE2 FTypeName log (const FTypeName& a) { FTypeName c(Op::myLog(a.val())); if (!a.depend()) return c; c.setDepend(a); const T& aval(a.val()); for(unsigned int i=0;i INLINE2 FTypeName log (const FTypeName& a) { FTypeName c(Op::myLog(a.val())); if (!a.depend()) return c; c.setDepend(a); const T& aval(a.val()); for(unsigned int i=0;i INLINE2 FTypeName sqrt (const FTypeName& a) { FTypeName c(Op::mySqrt(a.val())); if (!a.depend()) return c; T tmp(c.val()*Op::myTwo()); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName sqrt (const FTypeName& a) { FTypeName c(Op::mySqrt(a.val())); if (!a.depend()) return c; T tmp(c.val()*Op::myTwo()); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName sin (const FTypeName& a) { FTypeName c(Op::mySin(a.val())); if (!a.depend()) return c; T tmp(Op::myCos(a.val())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName sin (const FTypeName& a) { FTypeName c(Op::mySin(a.val())); if (!a.depend()) return c; T tmp(Op::myCos(a.val())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName cos (const FTypeName& a) { FTypeName c(Op::myCos(a.val())); if (!a.depend()) return c; T tmp(-Op::mySin(a.val())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName cos (const FTypeName& a) { FTypeName c(Op::myCos(a.val())); if (!a.depend()) return c; T tmp(-Op::mySin(a.val())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName tan (const FTypeName& a) { FTypeName c(Op::myTan(a.val())); if (!a.depend()) return c; T tmp(Op::myOne()+Op::mySqr(c.val())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName tan (const FTypeName& a) { FTypeName c(Op::myTan(a.val())); if (!a.depend()) return c; T tmp(Op::myOne()+Op::mySqr(c.val())); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName asin (const FTypeName& a) { FTypeName c(Op::myAsin(a.val())); if (!a.depend()) return c; T tmp(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val())))); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName asin (const FTypeName& a) { FTypeName c(Op::myAsin(a.val())); if (!a.depend()) return c; T tmp(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val())))); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName acos (const FTypeName& a) { FTypeName c(Op::myAcos(a.val())); if (!a.depend()) return c; T tmp(Op::myNeg(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val()))))); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName acos (const FTypeName& a) { FTypeName c(Op::myAcos(a.val())); if (!a.depend()) return c; T tmp(Op::myNeg(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val()))))); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName atan (const FTypeName& a) { FTypeName c(Op::myAtan(a.val())); if (!a.depend()) return c; T tmp(Op::myInv(Op::myOne()+Op::mySqr(a.val()))); c.setDepend(a); for(unsigned int i=0;i INLINE2 FTypeName atan (const FTypeName& a) { FTypeName c(Op::myAtan(a.val())); if (!a.depend()) return c; T tmp(Op::myInv(Op::myOne()+Op::mySqr(a.val()))); c.setDepend(a); for(unsigned int i=0;i struct Op< FTypeName > { typedef FTypeName T; typedef FTypeName Underlying; typedef typename Op::Base Base; static Base myInteger(const int i) { return Base(i); } static Base myZero() { return myInteger(0); } static Base myOne() { return myInteger(1);} static Base myTwo() { return myInteger(2); } static Base myPI() { return Op::myPI(); } static T myPos(const T& x) { return +x; } static T myNeg(const T& x) { return -x; } template static T& myCadd(T& x, const V& y) { return x+=y; } template static T& myCsub(T& x, const V& y) { return x-=y; } template static T& myCmul(T& x, const V& y) { return x*=y; } template static T& myCdiv(T& x, const V& y) { return x/=y; } static T myInv(const T& x) { return myOne()/x; } static T mySqr(const T& x) { return fadbad::sqr(x); } template static T myPow(const X& x, const Y& y) { return fadbad::pow(x,y); } static T mySqrt(const T& x) { return fadbad::sqrt(x); } static T myLog(const T& x) { return fadbad::log(x); } static T myExp(const T& x) { return fadbad::exp(x); } static T mySin(const T& x) { return fadbad::sin(x); } static T myCos(const T& x) { return fadbad::cos(x); } static T myTan(const T& x) { return fadbad::tan(x); } static T myAsin(const T& x) { return fadbad::asin(x); } static T myAcos(const T& x) { return fadbad::acos(x); } static T myAtan(const T& x) { return fadbad::atan(x); } static bool myEq(const T& x, const T& y) { return x==y; } static bool myNe(const T& x, const T& y) { return x!=y; } static bool myLt(const T& x, const T& y) { return xy; } static bool myGe(const T& x, const T& y) { return x>=y; } }; } // namespace fadbad #endif