C++11 Segmentation Fault with Boost Polynomials

I am trying to implement euclidian polynomials gcd using boost library. It compiles. However, I got segmentation fault on the runtime.

This is my code:

polynomial<double> gcd_euclid(polynomial<double> const &a, polynomial<double> const &b){

polynomial<double> s, old_t = zero_element(std::multiplies<polynomial<double>>());
polynomial<double> old_s, t = identity_element(std::multiplies<polynomial<double>>()); 
polynomial<double> r = b; 
polynomial<double> old_r = a;

while (r != zero_element(std::multiplies<polynomial<double>>())) {

    polynomial<double> quotient = old_r / r;

    old_r = r; 

    r = old_r - quotient * r;

    old_s = s;

    s = old_s - quotient * s;

    old_t = t;

    t = old_t - quotient * t;
}

return old_r;

}

Using gdb, it traces the segfault to this part in the polynomial.hpp:

 template <class U>
 polynomial& operator *=(const polynomial<U>& value)
 {
  // TODO: FIXME: use O(N log(N)) algorithm!!!
  polynomial const zero = zero_element(std::multiplies<polynomial>());
  if (value == zero)
  {
      *this = zero;
      return *this;
  }
  polynomial base(*this);
  this->multiplication(value[0]);
  for(size_type i = 1; i < value.size(); ++i)
  {
     polynomial t(base);
     t.multiplication(value[i]);
     size_type s = size() - i;
     for(size_type j = 0; j < s; ++j)
     {
        m_data[i+j] += t[j];
     }
     for(size_type j = s; j < t.size(); ++j)
        m_data.push_back(t[j]);
  }
  return *this;

}

specifically this line:

 m_data[i+j] += t[j];