diff --git a/src/bubble.cc b/src/bubble.cc
index 18b2018..59f4ee1 100644
--- a/src/bubble.cc
+++ b/src/bubble.cc
@@ -217,6 +217,73 @@ namespace ql
     res[2] = this->_czero;
   }
 
+    /*!
+   * Computes the Bubble derivative.
+   * Implementation of the formulae from Eq. (4.25) of A. Denner, arXiv:0709.1075.
+   *
+   * \param res output object res[0,1,2] the coefficients in the Laurent series
+   * \param mu2 is the squre of the scale mu
+   * \param m are the squares of the masses of the internal lines
+   * \param p are the four-momentum squared of the external lines
+   */
+  template<typename TOutput, typename TMass, typename TScale>
+  void Bubble<TOutput,TMass,TScale>::derivative(vector<TOutput> &res,
+                                                const TScale& mu2,
+                                                vector<TMass> const& m,
+                                                vector<TScale> const& p)
+  {
+    if (mu2 < 0)
+      throw RangeError("Bubble::derivative","mu2 is negative!");
+
+    if ((Real(m[0]) < 0 || Real(m[1]) < 0) || (Imag(m[0]) > 0 || Imag(m[1]) > 0))
+      throw RangeError("Bubble::derivative","Real masses must be positive, imag. negative");
+
+    if (res.size() != 3) { res.reserve(3); }
+    std::fill(res.begin(), res.end(), this->_czero);
+
+    if (this->iszero(p[0]))
+    {
+      if (this->iszero(Abs(m[0])) && this->iszero(Abs(m[1])))
+        return;
+      if (this->iszero(Abs(m[0])-Abs(m[1])))
+        res[0] = this->_cone / (6.0 * m[0]);
+    }
+    else
+    {
+      if (this->iszero(Abs(m[0])) && this->iszero(Abs(m[1])))
+        res[0] = - this->_cone / p[0];
+      else if (this->iszero(Min(m[0], m[1])))
+      {
+        TMass msq;
+        if (Abs(m[0]) >= Abs(m[1])) msq = m[0];
+        else msq = m[1];
+
+        if (this->iszero(Abs(p[0] - msq)))
+        {
+          res[1] = - this->_chalf / msq;
+          res[0] = - this->_chalf / msq * this->Lnrat(mu2, msq) - this->_cone / msq;
+        }
+        else
+          res[0] = - (this->_cone + msq / p[0] * this->Lnrat(msq - p[0], msq)) / p[0];
+      }
+      else
+      {
+        const TMass a = Sqrt(m[0] * m[1]);
+        const TMass c = a;
+        const TOutput b = m[0] + m[1] - TMass(p[0]) - this->_ieps;
+        const TOutput root = Sqrt(Pow(b, 2) - this->_cfour*a*c);
+        const TOutput sgn = TOutput(Sign(Real( Conjg(b) * root)));
+        const TOutput q = this->_chalf * (b + sgn*root);
+        const TOutput rm = q / a;
+        const TOutput r = rm;
+        res[0] = - this->_chalf * (m[0] - m[1]) / Pow(p[0], 2) * Log(m[1]/m[0]) +
+            Sqrt(m[0] * m[1]) / Pow(p[0], 2) * (this->_cone / r - r) * Log(r)
+            - (1.0 + (Pow(r, 2) + this->_cone) / (Pow(r, 2) - this->_cone) * Log(r)) / p[0];
+      }
+    }
+    return;
+  }
+
   // explicity tyoename declaration
   template class Bubble<complex,double,double>;
   template class Bubble<complex,complex,double>;
diff --git a/src/qcdloop/bubble.h b/src/qcdloop/bubble.h
index 1abf2bf..ace50db 100644
--- a/src/qcdloop/bubble.h
+++ b/src/qcdloop/bubble.h
@@ -43,5 +43,8 @@ namespace ql
 
     //! Special configuration I(0;m0,m1)
     void BB5(vector<TOutput> &res, TScale const& mu2, TMass const& m0, TMass const& m1) const;
+
+    //! Computes the Bubble derivative
+    void derivative(vector<TOutput> &res, TScale const& mu2, vector<TMass> const& m, vector<TScale> const& p);
     };
 }
diff --git a/src/qcdloop/wrapper.h b/src/qcdloop/wrapper.h
index 17cae90..593460a 100644
--- a/src/qcdloop/wrapper.h
+++ b/src/qcdloop/wrapper.h
@@ -55,6 +55,11 @@ extern"C" {
   qcomplex qli2q_(qdouble const& p1, qdouble const& m1, qdouble const& m2, qdouble const& mu2, int const& ep);
   qcomplex qli2qc_(qdouble const& p1, qcomplex const& m1, qcomplex const& m2, qdouble const& mu2, int const& ep);
 
+  complex  qli2p_(double const& p1, double const& m1, double const& m2, double const& mu2, int const& ep);
+  complex  qli2pc_(double const& p1, complex const& m1, complex const& m2, double const& mu2, int const& ep);
+  qcomplex qli2pq_(qdouble const& p1, qdouble const& m1, qdouble const& m2, qdouble const& mu2, int const& ep);
+  qcomplex qli2pqc_(qdouble const& p1, qcomplex const& m1, qcomplex const& m2, qdouble const& mu2, int const& ep);
+
   complex  qli3_(double const& p1, double const& p2, double const& p3, double const& m1, double const& m2, double const& m3, double const& mu2, int const& ep);
   complex  qli3c_(double const& p1, double const& p2, double const& p3, complex const& m1, complex const& m2, complex const& m3, double const& mu2, int const& ep);
   qcomplex qli3q_(qdouble const& p1, qdouble const& p2, qdouble const& p3, qdouble const& m1, qdouble const& m2, qdouble const& m3, qdouble const& mu2, int const& ep);
diff --git a/src/wrapper.cc b/src/wrapper.cc
index a167578..98febaa 100644
--- a/src/wrapper.cc
+++ b/src/wrapper.cc
@@ -533,7 +533,7 @@ void qlboxcq_(qcomplex (&out)[3], qdouble const& mu2, qcomplex const& m1, qcompl
     {
       std::cout << e.what() << std::endl;
       exit(-1);
-    }    
+    }
     return r[abs(ep)];
   }
 
@@ -588,6 +588,74 @@ void qlboxcq_(qcomplex (&out)[3], qdouble const& mu2, qcomplex const& m1, qcompl
     return rq[abs(ep)];
   }
 
+  complex qli2p_(double const& p1, double const& m1, double const& m2, double const& mu2, int const& ep)
+  {
+    try
+    {
+      mI2[0] = m1;
+      mI2[1] = m2;
+      pI2[0] = p1;
+      bb.derivative(r, mu2, mI2, pI2);
+    }
+    catch (std::exception &e)
+    {
+      std::cout << e.what() << std::endl;
+      exit(-1);
+    }
+    return r[abs(ep)];
+  }
+
+  complex qli2pc_(double const& p1, complex const& m1, complex const& m2, double const& mu2, int const& ep)
+  {
+    try
+    {
+      mI2c[0] = m1;
+      mI2c[1] = m2;
+      pI2[0] = p1;
+      bbc.derivative(r, mu2, mI2c, pI2);
+    }
+    catch (std::exception &e)
+    {
+      std::cout << e.what() << std::endl;
+      exit(-1);
+    }
+    return r[abs(ep)];
+  }
+
+  qcomplex qli2pq_(qdouble const& p1, qdouble const& m1, qdouble const& m2, qdouble const& mu2, int const& ep)
+  {
+    try
+    {
+      mI2q[0] = m1;
+      mI2q[1] = m2;
+      pI2q[0] = p1;
+      bbq.derivative(rq, mu2, mI2q, pI2q);
+    }
+    catch (std::exception &e)
+    {
+      std::cout << e.what() << std::endl;
+      exit(-1);
+    }
+    return rq[abs(ep)];
+  }
+
+  qcomplex qli2pqc_(qdouble const& p1, qcomplex const& m1, qcomplex const& m2, qdouble const& mu2, int const& ep)
+  {
+    try
+    {
+      mI2cq[0] = m1;
+      mI2cq[1] = m2;
+      pI2q[0] = p1;
+      bbcq.derivative(rq, mu2, mI2cq, pI2q);
+    }
+    catch (std::exception &e)
+    {
+      std::cout << e.what() << std::endl;
+      exit(-1);
+    }
+    return rq[abs(ep)];
+  }
+
   complex qli3_(double const& p1, double const& p2, double const& p3, double const& m1, double const& m2, double const& m3, double const& mu2, int const& ep)
   {
     try