// 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 _TADIFF_H #define _TADIFF_H #include #ifndef MaxLength #define MaxLength 40 #endif #include "fadbad.h" namespace fadbad { template class TValues { unsigned int m_n; U m_val[N]; public: TValues():m_n(0){std::fill(m_val,m_val+N,Op::myZero());} template explicit TValues(const V& val):m_n(1){m_val[0]=val;std::fill(m_val+1,m_val+N,Op::myZero());} U& operator[](const unsigned int i) { USER_ASSERT(i class TTypeNameHV // Heap Value { TValues m_val; mutable unsigned int m_rc; protected: virtual ~TTypeNameHV(){} public: TTypeNameHV():m_rc(0){} template explicit TTypeNameHV(const V& val):m_val(val),m_rc(0){} const U& val(const unsigned int i) const { return m_val[i]; } U& val(const unsigned int i) { return m_val[i]; } unsigned int length() const { return m_val.length(); } unsigned int& length() { return m_val.length(); } void decRef(TTypeNameHV*& pTTypeNameHV) const { if (--m_rc==0) { delete this; pTTypeNameHV=0;} } void incRef() const {++m_rc;} virtual void reset(){m_val.reset();} virtual unsigned int eval(const unsigned int k){return k+1;} }; template class TTypeName { private: struct SV // Stack Value refers to reference-counted Heap Value: { mutable TTypeNameHV* m_pTTypeNameHV; SV(TTypeNameHV* pTTypeNameHV):m_pTTypeNameHV(pTTypeNameHV){ m_pTTypeNameHV->incRef(); } SV(const typename TTypeName::SV& sv):m_pTTypeNameHV(sv.m_pTTypeNameHV){ m_pTTypeNameHV->incRef(); } ~SV(){ m_pTTypeNameHV->decRef(m_pTTypeNameHV); } TTypeNameHV* getTTypeNameHV() const { return m_pTTypeNameHV; } void setTTypeNameHV(TTypeNameHV* pTTypeNameHV) { if (m_pTTypeNameHV!=pTTypeNameHV) { m_pTTypeNameHV->decRef(m_pTTypeNameHV); m_pTTypeNameHV=pTTypeNameHV; m_pTTypeNameHV->incRef(); } } const U& val() const { return m_pTTypeNameHV->val(0); } const U& val(const unsigned int i) const { return m_pTTypeNameHV->val(i); } unsigned int length() const { return m_pTTypeNameHV->length(); } unsigned int& length() { return m_pTTypeNameHV->length(); } U& val(const unsigned int i) { return m_pTTypeNameHV->val(i); } void reset(){m_pTTypeNameHV->reset();} unsigned int eval(const unsigned int i){return m_pTTypeNameHV->eval(i);} } m_sv; public: typedef U UnderlyingType; TTypeName():m_sv(new TTypeNameHV()){} TTypeName(TTypeNameHV* pTTypeNameHV):m_sv(pTTypeNameHV){} explicit TTypeName(const typename TTypeName::SV& sv):m_sv(sv){} template /*explicit*/ TTypeName(const V& val):m_sv(new TTypeNameHV(val)){m_sv.length()=N;} TTypeName& operator=(const TTypeName& val) { if (this==&val) return *this; m_sv.setTTypeNameHV(val.m_sv.getTTypeNameHV()); return *this; } template TTypeName& operator=(const V& val) { m_sv.setTTypeNameHV(new TTypeNameHV(val)); m_sv.length()=N; return *this; } TTypeNameHV* getTTypeNameHV() const { return m_sv.getTTypeNameHV(); } void setTTypeNameHV(const TTypeNameHV* pTTypeNameHV) { m_sv.setTTypeNameHV(pTTypeNameHV); } const U& val() const { return m_sv.val(); } unsigned int length() const { return m_sv.length(); } const U& operator[](const unsigned int i) const { return m_sv.val(i); } U& operator[](const unsigned int i) { if (i>=m_sv.length()) m_sv.length()=i+1; return m_sv.val(i);} TTypeName& operator+=(const TTypeName& val); TTypeName& operator-=(const TTypeName& val); TTypeName& operator*=(const TTypeName& val); TTypeName& operator/=(const TTypeName& val); template TTypeName& operator+=(const V& val); template TTypeName& operator-=(const V& val); template TTypeName& operator*=(const V& val); template TTypeName& operator/=(const V& val); void reset(){m_sv.reset();} unsigned int eval(const unsigned int i){return m_sv.eval(i);} }; template bool operator==(const TTypeName& val1, const TTypeName& val2) { return Op::myEq(val1.val(),val2.val()); } template bool operator!=(const TTypeName& val1, const TTypeName& val2) { return Op::myNe(val1.val(),val2.val()); } template bool operator<(const TTypeName& val1, const TTypeName& val2) { return Op::myLt(val1.val(),val2.val()); } template bool operator<=(const TTypeName& val1, const TTypeName& val2) { return Op::myLe(val1.val(),val2.val()); } template bool operator>(const TTypeName& val1, const TTypeName& val2) { return Op::myGt(val1.val(),val2.val()); } template bool operator>=(const TTypeName& val1, const TTypeName& val2) { return Op::myGe(val1.val(),val2.val()); } template bool operator==(const TTypeName& val1, const V& val2) { return Op::myEq(val1.val(),val2); } template bool operator==(const V& val1, const TTypeName& val2) { return Op::myEq(val1,val2.val()); } template bool operator!=(const TTypeName& val1, const V& val2) { return Op::myNe(val1.val(),val2); } template bool operator!=(const V& val1, const TTypeName& val2) { return Op::myNe(val1,val2.val()); } template bool operator<(const TTypeName& val1, const V& val2) { return Op::myLt(val1.val(),val2); } template bool operator<(const V& val1, const TTypeName& val2) { return Op::myLt(val1,val2.val()); } template bool operator<=(const TTypeName& val1, const V& val2) { return Op::myLe(val1.val(),val2); } template bool operator<=(const V& val1, const TTypeName& val2) { return Op::myLe(val1,val2.val()); } template bool operator>(const TTypeName& val1, const V& val2) { return Op::myGt(val1.val(),val2); } template bool operator>(const V& val1, const TTypeName& val2) { return Op::myGt(val1,val2.val()); } template bool operator>=(const TTypeName& val1, const V& val2) { return Op::myGe(val1.val(),val2); } template bool operator>=(const V& val1, const TTypeName& val2) { return Op::myGe(val1,val2.val()); } // Binary operator base class: template class BinTTypeNameHV : public TTypeNameHV { TTypeNameHV* m_pOp1; TTypeNameHV* m_pOp2; public: BinTTypeNameHV(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):TTypeNameHV(val),m_pOp1(pOp1),m_pOp2(pOp2) { m_pOp1->incRef();m_pOp2->incRef(); } BinTTypeNameHV(TTypeNameHV* pOp1, TTypeNameHV* pOp2):TTypeNameHV(),m_pOp1(pOp1),m_pOp2(pOp2) { m_pOp1->incRef();m_pOp2->incRef(); } virtual ~BinTTypeNameHV() { m_pOp1->decRef(m_pOp1);m_pOp2->decRef(m_pOp2); } TTypeNameHV* op1() { return m_pOp1; } TTypeNameHV* op2() { return m_pOp2; } unsigned int op1Eval(const unsigned int k){return this->op1()->eval(k);} unsigned int op2Eval(const unsigned int k){return this->op2()->eval(k);} const U& op1Val(const unsigned int k) {return this->op1()->val(k);} const U& op2Val(const unsigned int k) {return this->op2()->val(k);} void reset(){op1()->reset();op2()->reset();TTypeNameHV::reset();} }; // Unary operator base class: template class UnTTypeNameHV : public TTypeNameHV { TTypeNameHV* m_pOp; public: UnTTypeNameHV(const U& val, TTypeNameHV* pOp):TTypeNameHV(val),m_pOp(pOp) { m_pOp->incRef(); } UnTTypeNameHV(TTypeNameHV* pOp):TTypeNameHV(),m_pOp(pOp) { m_pOp->incRef(); } virtual ~UnTTypeNameHV() { m_pOp->decRef(m_pOp); } TTypeNameHV* op() { return m_pOp; } unsigned int opEval(const unsigned int k){return this->op()->eval(k);} const U& opVal(const unsigned int k) {return this->op()->val(k);} void reset(){op()->reset();TTypeNameHV::reset();} }; // ADDITION: template struct TTypeNameADD : public BinTTypeNameHV { TTypeNameADD(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} TTypeNameADD(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); for(unsigned int i=this->length();ival(i)=this->op1Val(i)+this->op2Val(i); return this->length()=l; } private: void operator=(const TTypeNameADD&){} // not allowed }; template struct TTypeNameADD1 : public UnTTypeNameHV { const V m_a; TTypeNameADD1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} TTypeNameADD1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=m_a+this->opVal(0); this->length()=1; } for(unsigned int i=this->length();ival(i)=this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNameADD1&){} // not allowed }; template struct TTypeNameADD2 : public UnTTypeNameHV { const V m_b; TTypeNameADD2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} TTypeNameADD2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=this->opVal(0)+m_b; this->length()=1; } for(unsigned int i=this->length();ival(i)=this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNameADD2&){} // not allowed }; template TTypeName operator+(const TTypeName& val1, const TTypeName& val2) { TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? new TTypeNameADD(val1.val()+val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): new TTypeNameADD(val1.getTTypeNameHV(),val2.getTTypeNameHV()); return TTypeName(pHV); } /* template TTypeName operator+(const typename Op::Underlying& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameADD1::Underlying>(a+val2.val(), a, val2.getTTypeNameHV()): new TTypeNameADD1::Underlying>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator+(const TTypeName& val1, const typename Op::Underlying& b) { TTypeNameHV* pHV=val1.length()>0? new TTypeNameADD2::Underlying>(val1.val()+b, val1.getTTypeNameHV(), b): new TTypeNameADD2::Underlying>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } template TTypeName operator+(const typename Op::Base& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameADD1::Base>(a+val2.val(), a, val2.getTTypeNameHV()): new TTypeNameADD1::Base>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator+(const TTypeName& val1, const typename Op::Base& b) { TTypeNameHV* pHV=val1.length()>0? new TTypeNameADD2::Base>(val1.val()+b, val1.getTTypeNameHV(), b): new TTypeNameADD2::Base>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } */ template TTypeName operator+(const V& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameADD1(a+val2.val(), a, val2.getTTypeNameHV()): new TTypeNameADD1(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator+(const TTypeName& val1, const V& b) { TTypeNameHV* pHV=val1.length()>0? new TTypeNameADD2(val1.val()+b, val1.getTTypeNameHV(), b): new TTypeNameADD2(val1.getTTypeNameHV(), b); return TTypeName(pHV); } // SUBTRACTION: template struct TTypeNameSUB : public BinTTypeNameHV { TTypeNameSUB(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} TTypeNameSUB(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); for(unsigned int i=this->length();ival(i)=this->op1Val(i)-this->op2Val(i); return this->length()=l; } private: void operator=(const TTypeNameSUB&){} // not allowed }; template struct TTypeNameSUB1 : public UnTTypeNameHV { const V m_a; TTypeNameSUB1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} TTypeNameSUB1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=m_a-this->opVal(0); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myNeg(this->opVal(i)); return this->length()=l; } private: void operator=(const TTypeNameSUB1&){} // not allowed }; template struct TTypeNameSUB2 : public UnTTypeNameHV { const V m_b; TTypeNameSUB2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} TTypeNameSUB2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=this->opVal(0)-m_b; this->length()=1; } for(unsigned int i=this->length();ival(i)=this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNameSUB2&){} // not allowed }; template TTypeName operator-(const TTypeName& val1, const TTypeName& val2) { TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? new TTypeNameSUB(val1.val()-val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): new TTypeNameSUB(val1.getTTypeNameHV(),val2.getTTypeNameHV()); return TTypeName(pHV); } /* template TTypeName operator-(const typename Op::Underlying& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameSUB1::Underlying>(a-val2.val(), a, val2.getTTypeNameHV()): new TTypeNameSUB1::Underlying>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator-(const TTypeName& val1, const typename Op::Underlying& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameSUB2::Underlying>(val1.val()-b, val1.getTTypeNameHV(), b): new TTypeNameSUB2::Underlying>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } template TTypeName operator-(const typename Op::Base& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameSUB1::Base>(a-val2.val(), a, val2.getTTypeNameHV()): new TTypeNameSUB1::Base>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator-(const TTypeName& val1, const typename Op::Base& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameSUB2::Base>(val1.val()-b, val1.getTTypeNameHV(), b): new TTypeNameSUB2::Base>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } */ template TTypeName operator-(const V& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameSUB1(a-val2.val(), a, val2.getTTypeNameHV()): new TTypeNameSUB1(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator-(const TTypeName& val1, const V& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameSUB2(val1.val()-b, val1.getTTypeNameHV(), b): new TTypeNameSUB2(val1.getTTypeNameHV(), b); return TTypeName(pHV); } // MULTIPLICATION: template struct TTypeNameMUL : public BinTTypeNameHV { TTypeNameMUL(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} TTypeNameMUL(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=0;j<=i;++j) Op::myCadd(this->val(i),this->op1Val(j)*this->op2Val(i-j)); } return this->length()=l; } private: void operator=(const TTypeNameMUL&){} // not allowed }; template struct TTypeNameMUL1 : public UnTTypeNameHV { const V m_a; TTypeNameMUL1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} TTypeNameMUL1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=m_a*this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNameMUL1&){} // not allowed }; template struct TTypeNameMUL2 : public UnTTypeNameHV { const V m_b; TTypeNameMUL2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} TTypeNameMUL2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=this->opVal(i)*m_b; return this->length()=l; } private: void operator=(const TTypeNameMUL2&){} // not allowed }; template TTypeName operator*(const TTypeName& val1, const TTypeName& val2) { TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? new TTypeNameMUL(val1.val()*val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): new TTypeNameMUL(val1.getTTypeNameHV(),val2.getTTypeNameHV()); return TTypeName(pHV); } /* template TTypeName operator*(const typename Op::Underlying& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameMUL1::Underlying>(a*val2.val(), a, val2.getTTypeNameHV()): new TTypeNameMUL1::Underlying>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator*(const TTypeName& val1, const typename Op::Underlying& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameMUL2::Underlying>(val1.val()*b, val1.getTTypeNameHV(), b): new TTypeNameMUL2::Underlying>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } template TTypeName operator*(const typename Op::Base& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameMUL1::Base>(a*val2.val(), a, val2.getTTypeNameHV()): new TTypeNameMUL1::Base>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator*(const TTypeName& val1, const typename Op::Base& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameMUL2::Base>(val1.val()*b, val1.getTTypeNameHV(), b): new TTypeNameMUL2::Base>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } */ template TTypeName operator*(const V& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameMUL1(a*val2.val(), a, val2.getTTypeNameHV()): new TTypeNameMUL1(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator*(const TTypeName& val1, const V& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameMUL2(val1.val()*b, val1.getTTypeNameHV(), b): new TTypeNameMUL2(val1.getTTypeNameHV(), b); return TTypeName(pHV); } // DIVISION: template struct TTypeNameDIV : public BinTTypeNameHV { TTypeNameDIV(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} TTypeNameDIV(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); for(unsigned int i=this->length();ival(i)=this->op1Val(i); for(unsigned int j=1;j<=i;++j) Op::myCsub(this->val(i),this->op2Val(j)*this->val(i-j)); Op::myCdiv(this->val(i), this->op2Val(0)); } return this->length()=l; } private: void operator=(const TTypeNameDIV&){} // not allowed }; template struct TTypeNameDIV1 : public UnTTypeNameHV { const V m_a; TTypeNameDIV1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} TTypeNameDIV1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=m_a/this->opVal(0); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=1;j<=i;++j) Op::myCsub(this->val(i),this->opVal(j)*this->val(i-j)); Op::myCdiv(this->val(i), this->opVal(0)); } return this->length()=l; } private: void operator=(const TTypeNameDIV1&){} // not allowed }; template struct TTypeNameDIV2 : public UnTTypeNameHV { const V m_b; TTypeNameDIV2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} TTypeNameDIV2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=this->opVal(i)/m_b; return this->length()=l; } private: void operator=(const TTypeNameDIV2&){} // not allowed }; template TTypeName operator/(const TTypeName& val1, const TTypeName& val2) { TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? new TTypeNameDIV(val1.val()/val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): new TTypeNameDIV(val1.getTTypeNameHV(),val2.getTTypeNameHV()); return TTypeName(pHV); } /* template TTypeName operator/(const typename Op::Underlying& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameDIV1::Underlying>(a/val2.val(), a, val2.getTTypeNameHV()): new TTypeNameDIV1::Underlying>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator/(const TTypeName& val1, const typename Op::Underlying& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameDIV2::Underlying>(val1.val()/b, val1.getTTypeNameHV(), b): new TTypeNameDIV2::Underlying>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } template TTypeName operator/(const typename Op::Base& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameDIV1::Base>(a/val2.val(), a, val2.getTTypeNameHV()): new TTypeNameDIV1::Base>(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator/(const TTypeName& val1, const typename Op::Base& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameDIV2::Base>(val1.val()/b, val1.getTTypeNameHV(), b): new TTypeNameDIV2::Base>(val1.getTTypeNameHV(), b); return TTypeName(pHV); } */ template TTypeName operator/(const V& a, const TTypeName& val2) { TTypeNameHV* pHV=val2.length()>0 ? new TTypeNameDIV1(a/val2.val(), a, val2.getTTypeNameHV()): new TTypeNameDIV1(a, val2.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName operator/(const TTypeName& val1, const V& b) { TTypeNameHV* pHV=val1.length()>0 ? new TTypeNameDIV2(val1.val()/b, val1.getTTypeNameHV(), b): new TTypeNameDIV2(val1.getTTypeNameHV(), b); return TTypeName(pHV); } // COMPOUND ASSIGNMENTS: template TTypeName& TTypeName::operator+=(const TTypeName& val) { return (*this)=(*this)+val; } template TTypeName& TTypeName::operator-=(const TTypeName& val) { return (*this)=(*this)-val; } template TTypeName& TTypeName::operator*=(const TTypeName& val) { return (*this)=(*this)*val; } template TTypeName& TTypeName::operator/=(const TTypeName& val) { return (*this)=(*this)/val; } template template TTypeName& TTypeName::operator+=(const V& val) { return (*this)=(*this)+val; } template template TTypeName& TTypeName::operator-=(const V& val) { return (*this)=(*this)-val; } template template TTypeName& TTypeName::operator*=(const V& val) { return (*this)=(*this)*val; } template template TTypeName& TTypeName::operator/=(const V& val) { return (*this)=(*this)/val; } // UNARY MINUS template struct TTypeNameUMINUS : public UnTTypeNameHV { TTypeNameUMINUS(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} TTypeNameUMINUS(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=Op::myNeg(this->opVal(i)); return this->length()=l; } private: void operator=(const TTypeNameUMINUS&){} // not allowed }; template TTypeName operator-(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameUMINUS(Op::myNeg(val.val()),val.getTTypeNameHV()) : new TTypeNameUMINUS(val.getTTypeNameHV()); return TTypeName(pHV); } // UNARY PLUS template struct TTypeNameUPLUS : public UnTTypeNameHV { TTypeNameUPLUS(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} TTypeNameUPLUS(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=+this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNameUPLUS&){} // not allowed }; template TTypeName operator+(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameUPLUS(+val.val(),val.getTTypeNameHV()) : new TTypeNameUPLUS(val.getTTypeNameHV()); return TTypeName(pHV); } // POWER template struct TTypeNamePOW : public UnTTypeNameHV { TTypeNamePOW(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} TTypeNamePOW(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNamePOW&){} // not allowed }; template struct TTypeNamePOW1 : public UnTTypeNameHV { TTypeNamePOW1(const U& val, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2){} TTypeNamePOW1(TTypeNameHV* pOp2):UnTTypeNameHV(pOp2){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNamePOW1&){} // not allowed }; template struct TTypeNamePOW2 : public UnTTypeNameHV { TTypeNamePOW2(const U& val, TTypeNameHV* pOp1):UnTTypeNameHV(val,pOp1){} TTypeNamePOW2(TTypeNameHV* pOp1):UnTTypeNameHV(pOp1){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); for(unsigned int i=this->length();ival(i)=this->opVal(i); return this->length()=l; } private: void operator=(const TTypeNamePOW2&){} // not allowed }; template TTypeName pow(const TTypeName& val1, const TTypeName& val2) { TTypeName tmp(exp(val2*log(val1))); TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? new TTypeNamePOW(Op::myPow(val1.val(),val2.val()),tmp.getTTypeNameHV()) : new TTypeNamePOW(tmp.getTTypeNameHV()); return TTypeName(pHV); } /* template TTypeName pow(const typename Op::Underlying& a, const TTypeName& val2) { TTypeName tmp(exp(val2*Op::myLog(a))); TTypeNameHV* pHV=val2.length()>0 ? new TTypeNamePOW1::Underlying>(Op::myPow(a,val2.val()), tmp.getTTypeNameHV()) : new TTypeNamePOW1::Underlying>(tmp.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName pow(const TTypeName& val1, const typename Op::Underlying& b) { TTypeName tmp(exp(b*log(val1))); TTypeNameHV* pHV=val1.length()>0 ? new TTypeNamePOW2::Underlying>(Op::myPow(val1.val(),b), tmp.getTTypeNameHV()) : new TTypeNamePOW2::Underlying>(tmp.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName pow(const typename Op::Base& a, const TTypeName& val2) { TTypeName tmp(exp(val2*Op::Base>::myLog(a))); TTypeNameHV* pHV=val2.length()>0 ? new TTypeNamePOW1::Base>(Op::myPow(a,val2.val()), tmp.getTTypeNameHV()) : new TTypeNamePOW1::Base>(tmp.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName pow(const TTypeName& val1, const typename Op::Base& b) { TTypeName tmp(exp(b*log(val1))); TTypeNameHV* pHV=val1.length()>0 ? new TTypeNamePOW2::Base>(Op::myPow(val1.val(),b), tmp.getTTypeNameHV()) : new TTypeNamePOW2::Base>(tmp.getTTypeNameHV()); return TTypeName(pHV); } */ template TTypeName pow(const V& a, const TTypeName& val2) { TTypeName tmp(exp(val2*Op::myLog(a))); TTypeNameHV* pHV=val2.length()>0 ? new TTypeNamePOW1(Op::myPow(a,val2.val()), tmp.getTTypeNameHV()) : new TTypeNamePOW1(tmp.getTTypeNameHV()); return TTypeName(pHV); } template TTypeName pow(const TTypeName& val1, const V& b) { TTypeName tmp(exp(b*log(val1))); TTypeNameHV* pHV=val1.length()>0 ? new TTypeNamePOW2(Op::myPow(val1.val(),b), tmp.getTTypeNameHV()) : new TTypeNamePOW2(tmp.getTTypeNameHV()); return TTypeName(pHV); } // SQR template struct TTypeNameSQR : public UnTTypeNameHV { TTypeNameSQR(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} TTypeNameSQR(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=Op::mySqr(this->opVal(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); unsigned int m=(i+1)/2; for(unsigned int j=0;j::myCadd(this->val(i), this->opVal(i-j)*this->opVal(j)); Op::myCmul(this->val(i), Op::myTwo()); if (0==i%2) Op::myCadd(this->val(i), Op::mySqr(this->opVal(m))); } return this->length()=l; } private: void operator=(const TTypeNameSQR&){} // not allowed }; template TTypeName sqr(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameSQR(Op::mySqr(val.val()), val.getTTypeNameHV()) : new TTypeNameSQR(val.getTTypeNameHV()); return TTypeName(pHV); } // SQRT template struct TTypeNameSQRT : public UnTTypeNameHV { TTypeNameSQRT(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} TTypeNameSQRT(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=Op::mySqrt(this->opVal(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); unsigned int m=(i+1)/2; for(unsigned int j=1;j::myCadd(this->val(i), this->val(i-j)*this->val(j)); Op::myCmul(this->val(i), Op::myTwo()); if (0==i%2) Op::myCadd(this->val(i), Op::mySqr(this->val(m))); this->val(i)=(this->opVal(i)-this->val(i))/(Op::myTwo()*this->val(0)); } return this->length()=l; } private: void operator=(const TTypeNameSQRT&){} // not allowed }; template TTypeName sqrt(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameSQRT(Op::mySqrt(val.val()), val.getTTypeNameHV()): new TTypeNameSQRT(val.getTTypeNameHV()); return TTypeName(pHV); } // EXP template struct TTypeNameEXP : public UnTTypeNameHV { TTypeNameEXP(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} TTypeNameEXP(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=Op::myExp(this->opVal(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=0;j::myCadd(this->val(i), (Op::myOne()-Op::myInteger(j) / Op::myInteger(i))*this->opVal(i-j)*this->val(j)); } return this->length()=l; } private: void operator=(const TTypeNameEXP&){} // not allowed }; template TTypeName exp(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameEXP(Op::myExp(val.val()), val.getTTypeNameHV()): new TTypeNameEXP(val.getTTypeNameHV()); return TTypeName(pHV); } // LOG template struct TTypeNameLOG : public UnTTypeNameHV { TTypeNameLOG(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} TTypeNameLOG(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=Op::myLog(this->opVal(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=this->opVal(i); for(unsigned int j=1;j::myCsub(this->val(i), (Op::myOne()-Op::myInteger(j) / Op::myInteger(i))*this->opVal(j)*this->val(i-j)); Op::myCdiv(this->val(i), this->opVal(0)); } return this->length()=l; } private: void operator=(const TTypeNameLOG&){} // not allowed }; template TTypeName log(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameLOG(Op::myLog(val.val()), val.getTTypeNameHV()): new TTypeNameLOG(val.getTTypeNameHV()); return TTypeName(pHV); } // SIN template struct TTypeNameSIN : public UnTTypeNameHV { U m_COS[N]; TTypeNameSIN(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){m_COS[0]=Op::myCos(this->opVal(0));} TTypeNameSIN(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=Op::mySin(this->opVal(0)); m_COS[0]=Op::myCos(this->opVal(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=0;j::myCadd(this->val(i), Op::myInteger(j+1)*m_COS[i-1-j]*this->opVal(j+1)); Op::myCdiv(this->val(i), Op::myInteger(i)); m_COS[i]=Op::myZero(); for(unsigned int j=0;j::myCsub(m_COS[i], Op::myInteger(j+1)*this->val(i-1-j)*this->opVal(j+1)); Op::myCdiv(m_COS[i], Op::myInteger(i)); } return this->length()=l; } private: void operator=(const TTypeNameSIN&){} // not allowed }; template TTypeName sin(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameSIN(Op::mySin(val.val()), val.getTTypeNameHV()): new TTypeNameSIN(val.getTTypeNameHV()); return TTypeName(pHV); } // COS template struct TTypeNameCOS : public UnTTypeNameHV { U m_SIN[N]; TTypeNameCOS(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){m_SIN[0]=Op::mySin(this->opVal(0));} TTypeNameCOS(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} unsigned int eval(const unsigned int k) { unsigned int l=this->opEval(k); if (0==this->length()) { this->val(0)=Op::myCos(this->opVal(0)); m_SIN[0]=Op::mySin(this->opVal(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=0;j::myCsub(this->val(i), Op::myInteger(j+1)*m_SIN[i-1-j]*this->opVal(j+1)); Op::myCdiv(this->val(i), Op::myInteger(i)); m_SIN[i]=Op::myZero(); for(unsigned int j=0;j::myCadd(m_SIN[i], Op::myInteger(j+1)*this->val(i-1-j)*this->opVal(j+1)); Op::myCdiv(m_SIN[i], Op::myInteger(i)); } return this->length()=l; } private: void operator=(const TTypeNameCOS&){} // not allowed }; template TTypeName cos(const TTypeName& val) { TTypeNameHV* pHV=val.length()>0 ? new TTypeNameCOS(Op::myCos(val.val()), val.getTTypeNameHV()): new TTypeNameCOS(val.getTTypeNameHV()); return TTypeName(pHV); } // TAN template struct TTypeNameTAN : public BinTTypeNameHV { TTypeNameTAN(const U& val, TTypeNameHV* pOp, TTypeNameHV* pSqrCos):BinTTypeNameHV(val,pOp,pSqrCos){} TTypeNameTAN(TTypeNameHV* pOp, TTypeNameHV* pSqrCos):BinTTypeNameHV(pOp,pSqrCos){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); if (0==this->length()) { this->val(0)=Op::myTan(this->op1Val(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); this->val(i)=(this->op1Val(i)-this->val(i)/Op::myInteger(i))/this->op2Val(0); } return this->length()=l; } private: void operator=(const TTypeNameTAN&){} // not allowed }; template TTypeName tan(const TTypeName& val) { TTypeName tmp(sqr(cos(val))); TTypeNameHV* pHV=val.length()>0 ? new TTypeNameTAN(Op::myTan(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): new TTypeNameTAN(val.getTTypeNameHV(), tmp.getTTypeNameHV()); return TTypeName(pHV); } // ASIN template struct TTypeNameASIN : public BinTTypeNameHV { TTypeNameASIN(const U& val, TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(val,pOp,pSqrt){} TTypeNameASIN(TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(pOp,pSqrt){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); if (0==this->length()) { this->val(0)=Op::myAsin(this->op1Val(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); this->val(i)=(this->op1Val(i)-this->val(i)/Op::myInteger(i))/this->op2Val(0); } return this->length()=l; } private: void operator=(const TTypeNameASIN&){} // not allowed }; template TTypeName asin(const TTypeName& val) { TTypeName tmp(sqrt(Op::myOne()-sqr(val))); TTypeNameHV* pHV=val.length()>0 ? new TTypeNameASIN(Op::myAsin(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): new TTypeNameASIN(val.getTTypeNameHV(), tmp.getTTypeNameHV()); return TTypeName(pHV); } // ACOS template struct TTypeNameACOS : public BinTTypeNameHV { TTypeNameACOS(const U& val, TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(val,pOp,pSqrt){} TTypeNameACOS(TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(pOp,pSqrt){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); if (0==this->length()) { this->val(0)=Op::myAcos(this->op1Val(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); this->val(i)=Op::myNeg((this->op1Val(i)+this->val(i)/Op::myInteger(i))/this->op2Val(0)); } return this->length()=l; } private: void operator=(const TTypeNameACOS&){} // not allowed }; template TTypeName acos(const TTypeName& val) { TTypeName tmp(sqrt(Op::myOne()-sqr(val))); TTypeNameHV* pHV=val.length()>0 ? new TTypeNameACOS(Op::myAcos(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): new TTypeNameACOS(val.getTTypeNameHV(), tmp.getTTypeNameHV()); return TTypeName(pHV); } // ATAN template struct TTypeNameATAN : public BinTTypeNameHV { TTypeNameATAN(const U& val, TTypeNameHV* pOp, TTypeNameHV* p1pSqr):BinTTypeNameHV(val,pOp,p1pSqr){} TTypeNameATAN(TTypeNameHV* pOp, TTypeNameHV* p1pSqr):BinTTypeNameHV(pOp,p1pSqr){} unsigned int eval(const unsigned int k) { unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); if (0==this->length()) { this->val(0)=Op::myAtan(this->op1Val(0)); this->length()=1; } for(unsigned int i=this->length();ival(i)=Op::myZero(); for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); this->val(i)=(this->op1Val(i)-this->val(i)/Op::myInteger(i))/this->op2Val(0); } return this->length()=l; } private: void operator=(const TTypeNameATAN&){} // not allowed }; template TTypeName atan(const TTypeName& val) { TTypeName tmp(Op::myOne()+sqr(val)); TTypeNameHV* pHV=val.length()>0 ? new TTypeNameATAN(Op::myAtan(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): new TTypeNameATAN(val.getTTypeNameHV(), tmp.getTTypeNameHV()); return TTypeName(pHV); } // Ned's diff operator template struct DIFF : public UnTTypeNameHV { int m_b; DIFF(const U& val, TTypeNameHV* pOp, const int b):UnTTypeNameHV(val,pOp),m_b(b){} DIFF(TTypeNameHV* pOp, const int b):UnTTypeNameHV(pOp),m_b(b){} unsigned int eval(const unsigned int k) { // IN ORDER TO COMPUTE i'th ORDER COEFFICIENTS OF diff(m_o1,b) // WE NEED (i+b)'th ORDER COEFFICIENTS OF TaylorOp. unsigned int l=this->opEval(k+m_b); // NOW WE SHOULD HAVE op1 EXPANDED TO DEGREE (i+b)'th ORDER, // WE CAN PROCEED TO COMPUTE UP TO i'TH ORDER COEFFICIENTS // OF diff(op1,b). if (this->length()+m_blength();ii;--j){ fact*=j; } this->val(i)=this->opVal(i+m_b)*fact; } this->length()=l-m_b; } return this->length(); } private: void operator=(const DIFF&){} // not allowed }; template TTypeName diff(const TTypeName& val, const int b) { // IF THE ARGUMENT TAYLOR EXPANSION TaylorOp HAS BEEN EVALUATED TO // DEGREE i THEN WE CAN EVALUATE THE ZERO ORDER VALUE OF // diff(TaylorOp,i). // THIS FUNCTION EVALUATES THE 0.ORDER COEFFICIENT TTypeNameHV* pHV=0; if (val.length()>b) { unsigned int fact=1; for(unsigned int j=b;j>1;--j){ fact*=j; } pHV=new DIFF(val[b]*fact, val.getTTypeNameHV(), b); } else { pHV=new DIFF(val.getTTypeNameHV(), b); } return TTypeName(pHV); } template struct Op< TTypeName > { typedef TTypeName V; typedef TTypeName 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 V myPos(const V& x) { return +x; } static V myNeg(const V& x) { return -x; } template static V& myCadd(V& x, const Y& y) { return x+=y; } template static V& myCsub(V& x, const Y& y) { return x-=y; } template static V& myCmul(V& x, const Y& y) { return x*=y; } template static V& myCdiv(V& x, const Y& y) { return x/=y; } static V myInv(const V& x) { return myOne()/x; } static V mySqr(const V& x) { return fadbad::sqr(x); } template static V myPow(const X& x, const Y& y) { return fadbad::pow(x,y); } static V mySqrt(const V& x) { return fadbad::sqrt(x); } static V myLog(const V& x) { return fadbad::log(x); } static V myExp(const V& x) { return fadbad::exp(x); } static V mySin(const V& x) { return fadbad::sin(x); } static V myCos(const V& x) { return fadbad::cos(x); } static V myTan(const V& x) { return fadbad::tan(x); } static V myAsin(const V& x) { return fadbad::asin(x); } static V myAcos(const V& x) { return fadbad::acos(x); } static V myAtan(const V& x) { return fadbad::atan(x); } static bool myEq(const V& x, const V& y) { return x==y; } static bool myNe(const V& x, const V& y) { return x!=y; } static bool myLt(const V& x, const V& y) { return xy; } static bool myGe(const V& x, const V& y) { return x>=y; } }; } // namespace fadbad #endif