// GenFem - A high-level finite element library
// Copyright (C) 2010-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#ifndef _DFLOAT_H_
#define _DFLOAT_H_
#include <cmath>
#include <ostream>

template<class scalar=double> class dfloat
{

public:
  typedef scalar FPType;
  scalar v;
  scalar dv;
  dfloat(): v(scalar()),dv(scalar()) {}
  dfloat(const scalar v_): v(v_),dv(scalar()) {}
  dfloat(const scalar v_,const scalar dv_): v(v_),dv(dv_) {}
  dfloat(const dfloat<scalar>& df_): v(df_.v),dv(df_.dv) {}
  inline operator scalar() const {return v;} 
  inline const scalar& d() const {return dv;}
  inline const dfloat<scalar> operator+ (void) const {return dfloat<scalar>(v,dv);} //unary operator
  inline const dfloat<scalar> operator- (void) const {return dfloat<scalar>(-v,-dv);} //unary operator
  inline const dfloat<scalar> & operator+=(const dfloat<scalar> in) { v+=in.v;dv+=in.dv; return *this;}
  inline const dfloat<scalar> & operator-=(const dfloat<scalar> in) { v-=in.v;dv-=in.dv; return *this;}
  inline const dfloat<scalar> & operator*=(const dfloat<scalar> in) { dv*=in.v;dv+=v*in.dv;v*=in.v; return *this;}
  inline const dfloat<scalar> & operator/=(const dfloat<scalar> in) { dv*=in.v;dv-=v*in.dv;dv/=(in.v*in.v); v/=in.v; return *this;}
  inline const dfloat<scalar> & operator+=(const scalar in) { v+=in; return *this;}
  inline const dfloat<scalar> & operator-=(const scalar in) { v-=in; return *this;}
  inline const dfloat<scalar> & operator*=(const scalar in) { dv*=in;v*=in; return *this;}
  inline const dfloat<scalar> & operator/=(const scalar in) { dv/=in;v/=in; return *this;}
  inline const dfloat<scalar> operator+ (const dfloat<scalar> in) const {return dfloat<scalar>(*this)+=in;}
  inline const dfloat<scalar> operator- (const dfloat<scalar> in) const {return dfloat<scalar>(*this)-=in;}
  inline const dfloat<scalar> operator* (const dfloat<scalar> in) const {return dfloat<scalar>(*this)*=in;}
  inline const dfloat<scalar> operator/ (const dfloat<scalar> in) const {return dfloat<scalar>(*this)/=in;}
  inline const dfloat<scalar> operator+ (const scalar in) const {return dfloat<scalar>(*this)+=in;}
  inline const dfloat<scalar> operator- (const scalar in) const {return dfloat<scalar>(*this)-=in;}
  inline const dfloat<scalar> operator* (const scalar in) const {return dfloat<scalar>(*this)*=in;}
  inline const dfloat<scalar> operator/ (const scalar in) const {return dfloat<scalar>(*this)/=in;}
};


template <class scalar> std::ostream & operator << (std::ostream & output,const dfloat<scalar> & d)
{
  output << d.v << "|" << d.dv ;return output;
}

template<class scalar> const dfloat<scalar> operator+(const scalar l,const dfloat<scalar> r)  {return dfloat<scalar>(l+r.v,r.dv);}
template<class scalar> const dfloat<scalar> operator-(const scalar l,const dfloat<scalar> r)  {return dfloat<scalar>(l-r.v,-r.dv);}
template<class scalar> const dfloat<scalar> operator*(const scalar l,const dfloat<scalar> r)  {return dfloat<scalar>(l*r.v,l*r.dv);}
template<class scalar> const dfloat<scalar> operator/(const scalar l,const dfloat<scalar> r)  {return dfloat<scalar>(l/r.v,(-l*r.dv)/(r.v*r.v));}

template<class scalar> const dfloat<scalar> sqrt(const dfloat<scalar> arg)  {scalar sq=sqrt(arg.v);return dfloat<scalar>(sq,arg.dv / (2.0*sq));}
template<class scalar> const dfloat<scalar> fabs(const dfloat<scalar> arg)  {return dfloat<scalar>(fabs(arg.v),(arg.v>=0.0?arg.dv:-arg.dv));}
template<class scalar> const dfloat<scalar> pow(const dfloat<scalar> arg, const double e) {scalar pw=pow(arg.v,e);return dfloat<scalar>(pw,e*arg.dv*pow(arg.v,e-1.0));}
template<class scalar> const dfloat<scalar> log(const dfloat<scalar> arg)  {scalar lg=log(arg.v);return dfloat<scalar>(lg,arg.dv / arg.v);}
template<class scalar> const dfloat<scalar> log10(const dfloat<scalar> arg)  {scalar lg=log10(arg.v);return dfloat<scalar>(lg,arg.dv / (arg.v*log(10.0)));}
template<class scalar> const dfloat<scalar> exp(const dfloat<scalar> arg)  {scalar ex=exp(arg.v);return dfloat<scalar>(ex,arg.dv *ex);}
template<class scalar> const dfloat<scalar> sin(const dfloat<scalar> arg)  {scalar si=sin(arg.v);return dfloat<scalar>(si,cos(arg.v)*arg.dv);}
template<class scalar> const dfloat<scalar> cos(const dfloat<scalar> arg)  {scalar cs=cos(arg.v);return dfloat<scalar>(cs,-sin(arg.v)*arg.dv);}
template<class scalar> const dfloat<scalar> tan(const dfloat<scalar> arg)  {scalar tg=tan(arg.v);scalar sc=1.0/cos(arg.v); return dfloat<scalar>(tg,sc*sc*arg.dv);}
template<class scalar> const dfloat<scalar> asin(const dfloat<scalar> arg)  {scalar asi=asin(arg.v); return dfloat<scalar>(asi,arg.dv/(sqrt(1.0-arg.v*arg.v)));}
template<class scalar> const dfloat<scalar> acos(const dfloat<scalar> arg)  {scalar acs=acos(arg.v); return dfloat<scalar>(acs,-arg.dv/(sqrt(1.0-arg.v*arg.v)));}
template<class scalar> const dfloat<scalar> atan(const dfloat<scalar> arg)  {scalar atg=atan(arg.v); return dfloat<scalar>(atg,-arg.dv/(1.0+arg.v*arg.v));}
template<class scalar> const dfloat<scalar> atan2(const dfloat<scalar> argy,const dfloat<scalar> argx)
{
  scalar atg2=atan2(argy.v,argx.v);scalar r2=argy.v*argy.v+argx.v*argx.v; 
  return dfloat<scalar>(atg2,(argy.dv*argx.v-argx.dv*argy.v)/r2);
}
template<class scalar> const dfloat<scalar> sinh(const dfloat<scalar> arg)  {scalar si=sinh(arg.v);return dfloat<scalar>(si,cosh(arg.v)*arg.dv);}
template<class scalar> const dfloat<scalar> cosh(const dfloat<scalar> arg)  {scalar cs=cosh(arg.v);return dfloat<scalar>(cs,sinh(arg.v)*arg.dv);}
template<class scalar> const dfloat<scalar> tanh(const dfloat<scalar> arg)  {scalar tg=tanh(arg.v);scalar sc=1.0/cosh(arg.v); return dfloat<scalar>(tg,sc*sc*arg.dv);}
template<class scalar> const dfloat<scalar> asinh(const dfloat<scalar> arg)  {scalar asi=asinh(arg.v); return dfloat<scalar>(asi,arg.dv/(sqrt(1.0+arg.v*arg.v)));}
template<class scalar> const dfloat<scalar> acosh(const dfloat<scalar> arg)  {scalar acs=acosh(arg.v); return dfloat<scalar>(acs,-arg.dv/(sqrt(arg.v*arg.v-1.0)));}
template<class scalar> const dfloat<scalar> atanh(const dfloat<scalar> arg)  {scalar atg=atanh(arg.v); return dfloat<scalar>(atg,-arg.dv/(1-arg.v*arg.v));}

template <class func,class scalar> void gradient(func& f,scalar x,scalar y,scalar z,scalar tab[3])
{
  tab[0]=f(dfloat<scalar>(x,1.0),y,z).dv;
  tab[1]=f(x,dfloat<scalar>(y,1.0),z).dv;
  tab[2]=f(x,y,dfloat<scalar>(z,1.0)).dv;
}

template <class func,class scalar> const scalar laplacian(func& f,scalar x,scalar y,scalar z)
{
  scalar tab[3],lap=0.0;
  gradient(f,x,y,z,tab);
  for (int i=0;i<3;++i) lap+=tab[i]*tab[i];
  return lap;
}

template <class func1,class func2,class func3,class scalar> const scalar divergence(func1& fx,func2& fy,func3& fz,scalar x,scalar y,scalar z)
{
  return fx(dfloat<scalar>(x,1.0),y,z).dv+fy(x,dfloat<scalar>(y,1.0),z).dv+fz(x,y,dfloat<scalar>(z,1.0)).dv;
}

template <class func1,class func2,class func3,class scalar> void curl(func1& fx,func2& fy,func3& fz,scalar x,scalar y,scalar z, scalar tab[3])
{
  tab[0]=fz(x,dfloat<scalar>(y,1.0),z).dv-fy(x,y,dfloat<scalar>(z,1.0)).dv;
  tab[1]=fx(x,y,dfloat<scalar>(z,1.0)).dv-fz(dfloat<scalar>(x,1.0),y,z).dv;
  tab[2]=fy(dfloat<scalar>(x,1.0),y,z).dv-fx(x,dfloat<scalar>(y,1.0),z).dv;
}

#endif // _DFLOAT_H_