#include <boost/math/quaternion.hpp>
#include <boost/operators.hpp>


using namespace boost::math;

template<class T>
class Wrapper
  : boost::field_operators<Wrapper<T> >
{
public:
  Wrapper(const T& w = T())
    : v_(w)
  { }

  Wrapper<T>& operator+=(const Wrapper<T>& w)
  {
    v_ += w.v_;
    return *this;
  }

  Wrapper<T> operator+() const
  {
    return Wrapper<T>(+v_);
  }

  Wrapper<T>& operator-=(const Wrapper<T>& w)
  {
    v_ -= w.v_;
    return *this;
  }

  Wrapper<T> operator-()
  {
    return Wrapper<T>(-v_);
  }

  Wrapper<T>& operator*=(const Wrapper<T>& w)
  {
    v_ *= w.v_;
    return *this;
  }

  Wrapper<T>& operator/=(const Wrapper<T>& w)
  {
    v_ /= w.v_;
    return *this;
  }

private:
  T v_;
};

typedef Wrapper<int> wr;
typedef quaternion<wr> quat1;

typedef Wrapper<wr> wrr;
typedef quaternion<wrr> quat2;

int main(int argc, const char *argv[])
{
  // This works
  quat1 q1(wr(1), wr(2), wr(3), wr(4));
  pow(q1, 4);

  // This does not
  quat2 q2(wrr(3), wrr(4), wrr(5), wrr(6));
  //pow(q2, 4);

  return 0;
}