N维欧氏空间计算模板

Question

我一直在使用使用C++11存储坐标的std::vector[]代码。大多数情况下，这段代码使用2D或3D，但在我看来，它对于n维欧几里德距离的计算通常是有用的。

因为这些计算是在模拟中使用的，所以它们应该在不牺牲通用性的情况下尽可能快地实现。

// euclid.h
#include <vector>
#include <algorithm>
#include <stdexcept>

template <typename T>
std::vector<T> operator+(const std::vector<T>& a, const std::vector<T>& b)
{
    if (a.size() != b.size())
        throw std::domain_error("vector addition must have equal length vectors");
    std::vector<T> result;
    result.reserve(a.size());  
    std::transform(a.begin(), a.end(), b.begin(), 
                   std::back_inserter(result), std::plus<T>());
    return result;
}

template <typename T>
std::vector<T> operator-(const std::vector<T>& a, const std::vector<T>& b)
{
    if (a.size() != b.size())
        throw std::domain_error("vector subtraction must have equal length vectors");
    std::vector<T> result;
    result.reserve(a.size());
    std::transform(a.begin(), a.end(), b.begin(), 
                   std::back_inserter(result), std::minus<T>());
    return result;
}

template <typename T>
T squared_distance(const std::vector<T>& a, const std::vector<T>& b)
{
    if (a.size() != b.size())
        throw std::domain_error("squared_distance requires equal length vectors");
    return std::inner_product(a.begin(), a.end(), b.begin(), T(0), 
        std::plus<T>(), [](T x,T y){return (y-x)*(y-x);});

}

这里有一些驱动程序代码来演示使用。

// points.cpp
#include <iostream>
#include <vector>
#include "euclid.h"

template <typename T>
std::ostream& operator<<(std::ostream& out, const std::vector<T> &v)
{
    out << "{ ";
    for (auto p : v)
        out << p << ' ';
    return out << "}";
}

#define say(x) std::cout << "" #x " = " << (x) << std::endl
int main()
{
    std::vector<double> origin{0, 0, 0}, a{3, 4, 5}, b{-1, -2, -3},g{0,0,0,0};
    say(origin);
    say(a);
    say(b);
    say(a+b);
    say(a-b);
    say(b-a);
    say(squared_distance(origin,a));
    say(squared_distance(origin,b));
}

程序的输出如下所示：

origin = { 0 0 0 }
a = { 3 4 5 }
b = { -1 -2 -3 }
a+b = { 2 2 2 }
a-b = { 4 6 8 }
b-a = { -4 -6 -8 }
squared_distance(origin,a) = 50
squared_distance(origin,b) = 14

我对改进或批评现有代码的一般想法感兴趣，但我也有一些具体的问题：

1. squared_distance计算能更有效吗？
2. 在+和-操作符中有减少代码重复的好方法吗？
3. 我用int，float，double和std::complex做过测试。还有什么有用的？
4. 是否有必要容纳unsigned类型？
5. 我是否应该对可能的数值溢出或下溢做任何处理？
6. 在这个文件中有一个#ifndef头保护有什么意义吗？

Answer 1

为了减少冗余，我可能会将输入大小相同的检查单独移到函数模板中，然后从其他三个模板中调用它：

template <typename T>
void check_size(std::vector<T> const &a, std::vector<T> const &b) {
    if (a.size() != b.size())
        throw std::domain_error("vector addition must have equal length vectors");
}

然后，我至少按值传递一个输入参数，使用它作为计算的目标，并返回它：

template <typename T>
std::vector<T> operator+(std::vector<T> a, const std::vector<T>& b)
{
    check_size(a, b);
    std::transform(a.begin(), a.end(), 
                   b.begin(), 
                   a.begin(), 
                   std::plus<T>());
    return a;
}

虽然这并不总是提高速度，但它经常需要进行测试/分析，以了解在这种情况下它的工作效果如何。

另一个需要考虑的可能性是使用std::valarray代替。它是C++标准的一个被遗忘的子级(可以这么说)，但它是专门为支持快速的数值计算而设计的。它已经定义了+和-运算符以及sum和*运算符的等效项，因此您的代码如下所示：

std::valarray<double> origin{ 0, 0, 0 }, a{ 3, 4, 5 }, b{ -1, -2, -3 }, g{ 0, 0, 0, 0 };
say(origin);
say(a);
say(b);
say(a + b);
say(a - b);
say(b - a);
say(((a - origin) * (a - origin)).sum());
say(((b - origin) * (b - origin)).sum());

由于std::valarray定义了我们使用的所有操作符，我们所需要的唯一其他代码是say宏和operator<<重载(对其进行了详细的修改，以接受std::valarray参数而不是std::vector)。

考虑到你更具体的问题：

具体问题：

• squared_distance计算能更有效吗？

std::valarray和/或std::array可能会有所帮助--但对于相对较少的维度，我不希望看到太大的差异。

• 在+和-运算符中有减少代码重复的好方法吗？

上面的代码试图至少在某种程度上回答这个问题。

• 我已经用int、float、double和std::complex进行了测试。还有什么有用的？
• 是否有必要容纳无符号类型？

对于此任务，我认为无符号类型没有多大意义。是的，距离总是正数，但是你可以很容易地从中间计算中得到一个负数，而且你不想像用无符号类型那样把它包装成一个很大的数字。

• 我是否应该对可能的数值溢出或下溢做任何处理？

很难说。这基本上取决于你要使用的代码。如果你需要距离，有一些计算方法可以避免较大的平方距离(即使是中间距离)，而且计算速度也相当快--但它们仍然比计算平方距离慢，在这种情况下，最终结果通常也是你所处理的单个最大值，所以你不能很好地处理溢出。

在浮点中，减法是导致精度损失的主要原因之一。至少在理论上，如果可能的话，你可能更喜欢避免它，但是为了找到一个平方距离，我也不知道很多选择。

• 在这个文件中有一个#ifndef头保护有什么意义吗？

如果将编译器两次包含在同一个翻译单元中，编译器就会(或者应该)给出错误，这样就不会有任何伤害。当然，您不太可能直接在同一个翻译单元中包含它两次，但是通过另外两个标头将其包括在内的情况并不少见。

我也想提到std::array，但当我写这篇文章的时候，我发现Morwenn已经写了一篇关于这个的文章。

Answer 2

关于复制，您可以有一个用于应用函数元素的单一函数：

template< typename T, typename U, typename BinaryFunc >
auto elementwise_apply(const vector<T> &x, const vector<U> &y, BinaryFunc func)
     -> vector<decltype(func(x[0], y[0]))>
{
    // apply `func` to all arguments and return a vector
}

，然后根据它声明其他二进制操作。

template< typename T>
vector<T> operator+(const vector<T> &x, const vector<T> &y)
{
    return elementwise_apply(x, y, plus<T>());
}

以此类推。

我也同意看看是否可以使用array而不是vector --在编译时(而不是在运行时)可以实现的越多，越好。同样，更少的动态分配和更少的非内联函数通常会导致更高效的代码。

我还建议您考虑创建一个具有point成员(或vector成员)的array结构，而不是使用裸露的array，这样您的所有函数都只适用于实际引用点的对象。

Answer 3

我知道这不会回答您的任何问题，但是您应该考虑使用std::vector<T>来表示数据，而不是使用std::array<T, N>。这就是为什么这个设计可能更有趣的原因：

• std::array<T, N>使用堆栈内存，而std::vector<T>使用堆内存。我敢打赌，堆栈分配的数组比堆分配的向量要高效一些。
• 从数学上讲，改变向量的大小是毫无意义的。因此，固定大小的容器也将确保大小不会改变。
• 你不用再为reserve操心了。
• 在函数中，要确保所加或减的向量具有相同的大小。如果使用std::array<T, N>存储数据，将得到隐式编译时检查，而不是显式运行时检查:模板 std::array operator+(const std::array& a、const std::array& b) { std::array结果；for (size_t i=0；i