diff --git a/lib/Function.cpp b/lib/Function.cpp
index 09c1873..a546d8d 100644
--- a/lib/Function.cpp
+++ b/lib/Function.cpp
@@ -100,6 +100,50 @@ double DoubleFunction::operator()(void) const
     return (*this)(nullptr);
 }
 
+// arithmetic operators ////////////////////////////////////////////////////////
+DoubleFunction DoubleFunction::operator-(void) const
+{
+    DoubleFunction resFunc;
+    
+    auto res = [this](const double *arg){return -(*this)(arg);};
+    resFunc.setFunction(res, getNArg());
+    
+    return resFunc;
+}
+
+#define MAKE_SELF_FUNC_OP(op)\
+DoubleFunction & DoubleFunction::operator op##=(const DoubleFunction &f)\
+{\
+    DoubleFunction copy(*this);\
+    checkSize(f.getNArg());\
+    auto res = [&f, copy](const double *arg){return copy(arg) op f(arg);};\
+    setFunction(res, getNArg());\
+    return *this;\
+}\
+DoubleFunction & DoubleFunction::operator op##=(const DoubleFunction &&f)\
+{\
+    *this op##= f;\
+    return *this;\
+}
+
+#define MAKE_SELF_SCALAR_OP(op)\
+DoubleFunction & DoubleFunction::operator op##=(const double x)\
+{\
+    DoubleFunction copy(*this);\
+    auto res = [x, copy](const double *arg){return copy(arg) op x;};\
+    setFunction(res, getNArg());\
+    return *this;\
+}\
+
+MAKE_SELF_FUNC_OP(+)
+MAKE_SELF_FUNC_OP(-)
+MAKE_SELF_FUNC_OP(*)
+MAKE_SELF_FUNC_OP(/)
+MAKE_SELF_SCALAR_OP(+)
+MAKE_SELF_SCALAR_OP(-)
+MAKE_SELF_SCALAR_OP(*)
+MAKE_SELF_SCALAR_OP(/)
+
 /******************************************************************************
  *                    DoubleFunctionSample implementation                     *
  ******************************************************************************/
diff --git a/lib/Function.hpp b/lib/Function.hpp
index dd37374..58b1f88 100644
--- a/lib/Function.hpp
+++ b/lib/Function.hpp
@@ -35,6 +35,7 @@ BEGIN_NAMESPACE
 class DoubleFunction
 {
 private:
+    // function type
     typedef std::function<double(const double *)> vecFunc;
 public:
     // constructor
@@ -53,6 +54,20 @@ public:
     double operator()(void) const;
     template <typename... Ts>
     double operator()(const double arg0, const Ts... args) const;
+    // arithmetic operators
+    DoubleFunction   operator-(void) const;
+    DoubleFunction & operator+=(const DoubleFunction &f);
+    DoubleFunction & operator+=(const DoubleFunction &&f);
+    DoubleFunction & operator-=(const DoubleFunction &f);
+    DoubleFunction & operator-=(const DoubleFunction &&f);
+    DoubleFunction & operator*=(const DoubleFunction &f);
+    DoubleFunction & operator*=(const DoubleFunction &&f);
+    DoubleFunction & operator/=(const DoubleFunction &f);
+    DoubleFunction & operator/=(const DoubleFunction &&f);
+    DoubleFunction & operator+=(const double x);
+    DoubleFunction & operator-=(const double x);
+    DoubleFunction & operator*=(const double x);
+    DoubleFunction & operator/=(const double x);
 private:
     // error checking
     void checkSize(const Index nPar) const;
@@ -61,6 +76,9 @@ private:
     vecFunc               f_;
 };
 
+/******************************************************************************
+ *                  DoubleFunction template implementation                    *
+ ******************************************************************************/
 template <typename... Ts>
 double DoubleFunction::operator()(const double arg0, const Ts... args) const
 {
@@ -74,6 +92,63 @@ double DoubleFunction::operator()(const double arg0, const Ts... args) const
     return (*this)(arg);
 }
 
+/******************************************************************************
+ *                DoubleFunction inline arithmetic operators                  *
+ ******************************************************************************/
+#define MAKE_INLINE_FUNC_OP(op)\
+inline DoubleFunction operator op(DoubleFunction lhs,\
+                                  const DoubleFunction &rhs)\
+{\
+    lhs op##= rhs;\
+    return lhs;\
+}\
+inline DoubleFunction operator op(DoubleFunction lhs,\
+                                  const DoubleFunction &&rhs)\
+{\
+    return lhs op rhs;\
+}
+
+#define MAKE_INLINE_RSCALAR_OP(op)\
+inline DoubleFunction operator op(DoubleFunction lhs, const double rhs)\
+{\
+    lhs op##= rhs;\
+    return lhs;\
+}\
+
+#define MAKE_INLINE_LSCALAR_OP(op)\
+inline DoubleFunction operator op(const double lhs, DoubleFunction rhs)\
+{\
+    rhs op##= lhs;\
+    return rhs;\
+}
+
+MAKE_INLINE_FUNC_OP(+)
+MAKE_INLINE_FUNC_OP(-)
+MAKE_INLINE_FUNC_OP(*)
+MAKE_INLINE_FUNC_OP(/)
+MAKE_INLINE_RSCALAR_OP(+)
+MAKE_INLINE_RSCALAR_OP(-)
+MAKE_INLINE_RSCALAR_OP(*)
+MAKE_INLINE_RSCALAR_OP(/)
+MAKE_INLINE_LSCALAR_OP(+)
+MAKE_INLINE_LSCALAR_OP(*)
+
+// special case for scalar - function
+inline DoubleFunction operator-(const double lhs, DoubleFunction rhs)
+{
+    return (-rhs) + lhs;
+}
+
+// special case for scalar/function
+inline DoubleFunction operator/(const double lhs, DoubleFunction rhs)
+{
+    auto res = [lhs, &rhs](const double *arg){return lhs/rhs(arg);};
+    
+    rhs.setFunction(res, rhs.getNArg());
+    
+    return rhs;
+}
+
 /******************************************************************************
  *                      DoubleFunctionSample class                            *
  ******************************************************************************/
@@ -95,7 +170,7 @@ public:
 
 template <typename... Ts>
 DSample DoubleFunctionSample::operator()(const double arg0,
-                                        const Ts... args) const
+                                         const Ts... args) const
 {
     const double arg[] = {arg0, args...};