我应该如何使用表达式模板来实现数学向量类的标量乘法

How should I use expression templates in order to implement scalar multiplication for a mathematical vector class

本文关键字:向量 标量 实现 何使用 表达式 我应该      更新时间:2023-10-16

请考虑数学向量类的以下(部分)实现(基本上是您可以在维基百科有关表达式模板的文章中找到的代码):

namespace math
{
    template<class E>
    class vector_expression
    {
    public:
        std::size_t size() const {
            return static_cast<E const&>(*this).size();
        }
        double operator[](size_t i) const
        {
            if (i >= size())
                throw std::length_error("");
            return static_cast<E const&>(*this)[i];
        }
        operator E&() { return static_cast<E&>(*this); }
        operator E const&() const { return static_cast<E const&>(*this); }
    }; // class vector_expression
    template<class E1, class E2>
    class vector_sum
        : public vector_expression<vector_sum<E1, E2>>
    {
    public:
        vector_sum(vector_expression<E1> const& e1, vector_expression<E2> const& e2)
            : m_e1(e1), m_e2(e2)
        {
            if (e1.size() != e2.size())
                throw std::logic_error("");
        }
        std::size_t size() const {
            return m_e1.size(); // == m_e2.size()
        }
        double operator[](std::size_t i) const {
            return m_e1[i] + m_e2[i];
        }
    private:
        E1 const& m_e1;
        E2 const& m_e2;
    }; // class vector_sum
    template<typename E1, typename E2>
    vector_sum<E1, E2> operator+(vector_expression<E1> const& e1, vector_expression<E2> const& e2) {
        return { e1, e2 };
    }
    template<typename T>
    class vector
        : public vector_expression<vector<T>>
    {
    public:
        vector(std::size_t d)
            : m_data(d)
        { }
        vector(std::initializer_list<T> init)
            : m_data(init)
        { }
        template<class E>
        vector(vector_expression<E> const& expression)
            : m_data(expression.size())
        {
            for (std::size_t i = 0; i < expression.size(); ++i)
                m_data[i] = expression[i];
        }
        std::size_t size() const {
            return m_data.size();
        }
        double  operator[](size_t i) const { return m_data[i]; }
        double& operator[](size_t i) { return m_data[i]; }
    private:
        std::vector<T> m_data;
    }; // class vector    
} // namespace math

我应该如何扩展此实现以允许以下操作:

vector<double> x = { ... };
auto y = 4711 * x; // or y = x * 4711
auto z = 1 + x; // or x + 1, which should yield z[i] = x[i] + 1

我想我需要类似的东西

namespace math
{
    template<class E, typename T>
    class vector_product
        : public vector_expression<vector_product<E, T>>
    {
    public:
        vector_product(vector_expression<E> const& e, T const& t)
            : m_e(e), m_t(t)
        { }
        std::size_t size() const {
            return m_e.size();
        }
        double operator[](std::size_t i) const {
            return m_e[i] * m_t;
        }
    private:
        E const& m_e;
        T const& m_t;
    }; // class vector_product
    template<class E, typename T>
    vector_product<E, T> operator*(vector_expression<E> const& e, T const& t) {
        return { e, t };
    }
    template<class E, typename T>
    vector_product<E, T> operator*(T const& t, vector_expression<E> const& e) {
        return e * t;
    }   
} // namespace math

但我不知道这是否是一个好方法。那么,我应该怎么做呢?我应该添加任何复制或移动构造函数/赋值运算符吗?我想不是,因为隐式的应该做得很好,因为vector的唯一成员变量是 STL 类型。

我只是简单地扩展vector_sum以允许E2 double,如果是的话,可以优雅地处理它。这将涉及在构造函数中获取E1 const&E2 const&的参数,可能不会抛出大小差异(因为标量没有大小),以及重写operator[]不执行索引。对于最后一部分,如下所示:

    double operator[](std::size_t i) const {
        return m_e1[i] + get(m_e2, i);
    }
private:
    template <class E>
    double get(E const& rhs, std::size_t i) const {
        return rhs[i];
    }
    double get(double scalar, std::size_t ) const {
        return scalar;
    }

这样,如果您要添加两个vector_expression,您将执行索引,但如果您要添加vector_expressiondouble - 您甚至不会尝试索引到double中。此切换发生在编译时,因此没有运行时开销。

然后,您只需要再添加几个operator+

template <typename E1>
vector_sum<E1, double> operator+(vector_expression<E1> const& e1, double d) {
        return {e1, d};
}
template <typename E1>
vector_sum<E1, double> operator+(double d, vector_expression<E1> const& e1) {
        return {e1, d};
}

这使您可以编写:

math::vector<int> x = {1, 2, 3, 4};
math::vector<int> y = {2, 3, 4, 5};    
auto sum = 3 + x + 1;
<小时 />

保留对 const 的引用可能不是正确的做法 - 如果您a+b+c这样做,您最终将保留对临时a+b的引用。您可能只想保留对实际vector的引用,并保留所有中间对象的副本。

<小时 />

为了支持vector_product,您可能希望vector_sum<E1,E2>真正vector_binary_op<E1,E2,std::plus<>>然后vector_product<E1,E2>应该vector_binary_op<E1,E2,std::multiplies<>>。这样,您就不会有所有的重复。

受到 Yakk 和 Barry 的启发,我终于想出了以下内容:

namespace math
{
    template<class E>
    class expression
    {
    public:
        auto size() const {
            return static_cast<E const&>(*this).size();
        }
        auto operator[](std::size_t i) const
        {
            if (i >= size())
                throw std::length_error("");
            return static_cast<E const&>(*this)[i];
        }
        operator E&() { return static_cast<E&>(*this); }
        operator E const&() const { return static_cast<E const&>(*this); }
    }; // class expression
    template<typename T, class Allocator = std::allocator<T>>
    class vector
        : public expression<vector<T>>
    {
    private:
        using data_type = std::vector<T, Allocator>;
        data_type m_data;
    public:
        using value_type = T;
        using allocator_type = Allocator;
        using size_type = typename data_type::size_type;
        using difference_type = typename data_type::difference_type;
        using reference = typename data_type::reference;
        using const_reference = typename data_type::const_reference;
        using pointer = typename data_type::pointer ;
        using const_pointer = typename data_type::const_pointer;
        vector(size_type d)
            : m_data(d)
        { }
        vector(std::initializer_list<value_type> init)
            : m_data(init)
        { }
        template<class E>
        vector(expression<E> const& expression)
            : m_data(expression.size())
        {
            for (size_type i = 0; i < expression.size(); ++i)
                m_data[i] = expression[i];
        }
        size_type size() const {
            return m_data.size();
        }
        value_type  operator[](size_type i) const { return m_data[i]; }
        value_type& operator[](size_type i)       { return m_data[i]; };
    }; // class vector
    namespace detail
    {
        template<typename T>
        class scalar
            : public expression<scalar<T>>
        {
        public:
            using value_type = T;
            using allocator_type = std::allocator<void>;
            using size_type = typename std::allocator<T>::size_type;
            using difference_type = typename std::allocator<T>::difference_type;
            using reference = typename std::allocator<T>::reference;
            using const_reference = typename std::allocator<T>::const_reference;
            using pointer = typename std::allocator<T>::pointer;
            using const_pointer = typename std::allocator<T>::const_pointer;
            scalar(value_type value)
                : m_value(value)
            { }
            size_type size() const {
                return 0;
            }
            operator value_type&() { return static_cast<value_type&>(*this); }
            operator value_type const&() const { return static_cast<value_type const&>(*this); }
            value_type  operator[](size_type i) const { return m_value; }
            value_type& operator[](size_type i) { return m_value; }
        private:
            value_type m_value;
        }; // class scalar
        template<class>
        struct is_scalar : std::false_type { };
        template<class T>
        struct is_scalar<scalar<T>> : std::true_type { };
    } // namespace detail
    template<class E1, class E2, class BinaryOperation>
    class vector_binary_operation
        : public expression<vector_binary_operation<E1, E2, BinaryOperation>>
    {
    public:
        using value_type = decltype(BinaryOperation()(typename E1::value_type(), typename E2::value_type()));
        using allocator_type = std::conditional_t<
            detail::is_scalar<E1>::value,
            typename E2::allocator_type::template rebind<value_type>::other,
            typename E1::allocator_type::template rebind<value_type>::other>;
    private:
        using vector_type = vector<value_type, allocator_type>;
    public:
        using size_type = typename vector_type::size_type;
        using difference_type = typename vector_type::difference_type;
        using reference = typename vector_type::reference;
        using const_reference = typename vector_type::const_reference;
        using pointer = typename vector_type::pointer;
        using const_pointer = typename vector_type::const_pointer;
        vector_binary_operation(expression<E1> const& e1, expression<E2> const& e2, BinaryOperation op)
            : m_e1(e1), m_e2(e2),
              m_op(op)
        {
            if (e1.size() > 0 && e2.size() > 0 && !(e1.size() == e2.size()))
                throw std::logic_error("");
        }
        size_type size() const {
            return m_e1.size(); // == m_e2.size()
        }
        value_type operator[](size_type i) const {
            return m_op(m_e1[i], m_e2[i]);
        }
    private:
        E1 m_e1;
        E2 m_e2;
        //E1 const& m_e1;
        //E2 const& m_e2;
        BinaryOperation m_op;
    }; // class vector_binary_operation
    template<class E1, class E2>
    vector_binary_operation<E1, E2, std::plus<>>
    operator+(expression<E1> const& e1, expression<E2> const& e2) {
        return{ e1, e2, std::plus<>() };
    }
    template<class E1, class E2>
    vector_binary_operation<E1, E2, std::minus<>>
    operator-(expression<E1> const& e1, expression<E2> const& e2) {
        return{ e1, e2, std::minus<>() };
    }
    template<class E1, class E2>
    vector_binary_operation<E1, E2, std::multiplies<>>
    operator*(expression<E1> const& e1, expression<E2> const& e2) {
        return{ e1, e2, std::multiplies<>() };
    }
    template<class E1, class E2>
    vector_binary_operation<E1, E2, std::divides<>>
    operator/(expression<E1> const& e1, expression<E2> const& e2) {
        return{ e1, e2, std::divides<>() };
    }
    template<class E, typename T>
    vector_binary_operation<E, detail::scalar<T>, std::divides<>>
    operator/(expression<E> const& expr, T val) {
        return{ expr, detail::scalar<T>(val), std::divides<>() };
    }
    template<class E, typename T>
    vector_binary_operation<E, detail::scalar<T>, std::multiplies<>>
    operator*(T val, expression<E> const& expr) {
        return{ expr, detail::scalar<T>(val), std::multiplies<>() };
    }
    template<class E, typename T>
    vector_binary_operation<E, detail::scalar<T>, std::multiplies<>>
    operator*(expression<E> const& expr, T val) {
        return{ expr, detail::scalar<T>(val), std::multiplies<>() };
    }
} // namespace math

这允许对同一size的两个vector的运算+, -, /, *,以及vector x 和值a的乘法a * xx * a。此外,我们可以划分x / a ,但不能a / x(因为这没有意义)。我认为这是最合理的解决方案。

但是,仍然存在一些问题:我已将Allocator模板参数添加到vector类。在vector_binary_operation我需要知道生成的vector类型。两个expression都有可能有不同的allocator_type。我决定选择vector_binary_operation第一expressionallocator_type。我不认为这是一个真正的问题,因为我认为在这种情况下使用不同的Allocator没有多大意义。

更大的问题是我不知道我需要如何处理vector_binary_operation中的expression成员变量。将它们声明为 const 引用是有意义的,因为代码的全部意义在于避免不必要的副本。但是,正如巴里指出的那样,如果我们写sum = a + b + c,我们最终会保留对临时a + b的引用。做sum[0]会暂时operator()[0]。但是该对象在上一行之后被删除。

我不知道我需要在这里做什么,并问了另一个问题。

E=double专门化template<class E> class vector_expression

添加最小/最大大小的概念(因为双精度是 0 到无限维度),也许is_scalar(可能不需要,除了投射到标量?

杀死vector_sum并拿走它的玩具。 制作对元素进行元素操作的binop_vector

friend binop_vector<vector_expression,R,ElemAdd>
operator+( vector_expression const& l. R const& r ){
  return {l,r};
}
friend binop_vector<vector_expression<double>,vector_expression,ElemAdd>
operator+( double const& l. vector_expression const const& r ){
  return {l,r};
}

对于*来说,ElemMult vector_expression类似. 这 binop 使用元素明智的操作来实现[]

更改断言以确保最小/最大大小重叠。 在binop_vector中报告交集。

vector_expression<double>重载的最小值为 0 最大值 -1 (size_t),并始终返回其值。 如果需要多个标量类型,请编写scalar_expression<T>并从中继承vector_expression<double>(等)。

以上没有经过测试,只是在我坐在车库里时写在手机上。