2 个向量 C++ 的点积

Question

我必须编写输出两个向量的点积的程序。

仅使用 Double 类型组织计算以获得尽可能准确的结果。

输入应如下所示：

   N - vector length
   x1, x2,..., xN      co-ordinates of vector x (double type)
   y1, y2,..., yN      co-ordinates of vector y (double type)

输入样本：

   4
   1.0e20 -1.0e3 0.1 1.0e20
   1.0 4.0 -4.0 -1.0

以上向量的输出：

   -4000.4 

还有我的代码（我还没有使用 cin 因为起初我想用示例输入编写工作程序）：

#include <iostream>
#include <numeric>
#include <vector>
#include <functional>

int main(){
//double N; //length of both vectors , will be used when I will have to input vectors by cin
//std::cin >> N;
//N = 4;

std::vector<double> x{1.0e20, -1.0e3, 0.1, 1.0e20};
std::vector<double> y{1.0, 4.0, -4.0, -1.0};
double result = std::inner_product(x.begin(), x.end(), y.begin(), 0);

std::cout << result;
return 0;
}

我的输出是 -2.14748e+09 所以它甚至不接近预期的输出。我应该怎么做才能使它工作？

Answer 1

（第一）问题

这是 中内积的函数模板<numeric>：

template <class InputIterator1, class InputIterator2, class T>
   T inner_product (InputIterator1 first1, InputIterator1 last1,
                    InputIterator2 first2, T init);

请注意，定义T输出类型的是init参数。因此，鉴于您的输入：

std::inner_product(x.begin(), x.end(), y.begin(), 0);

init = 0，因此类型T为int。因此，当算法运行时，它会将double值类型转换为ints，最终将返回一个未定义的int值。

“修复”和第二个问题

要解决这个问题，您所要做的就是给出一个正确类型的init值（即，给出 adouble作为init参数）。只要0.0会做：

std::inner_product(x.begin(), x.end(), y.begin(), 0.0);

现在，当您使用该修复程序编译并运行程序时，它仍然会输出错误的结果：0

这是因为当inner_product函数累加值时，它使用标准double加法。因此，您会受到标准double不精确性的影响，它的机器 epsilon为 2^(-52) — 2.22E-16 或小数点后第 16 位的不精确性 — 这意味着，对于数字 1E20，(1E20 + x) = 1E20 对于所有 x < 2^(-52)*1E20 ？22204.46。

为了说明这一点，让我们1E20 + 23000在python解释器中添加（提醒python使用IEEE-754浮点运算，它等于double标准C++编译器中的精度）：

>>> 1e20 + 23000
1.0000000000000002e+20

因此，您会看到添加中忽略/“吸收”了少于两万的任何内容。

由于您的其他数字小于 22204.46，因此 1e20 只会“吸收”它们，直到将其添加到 -1E20，然后“取消”并返回0.

（简单）修复

解决第二个问题的最简单方法是使用long double而不是double. 这种更精确的双精度类型的机器 epsilon 为 2^(-63) — 1.08E-19 或大约 19 个小数位 — 这意味着，对于您的输入 1E20，不精确将等于 2^(-63) *1E20，或约 10.84。运行该程序，输出将为-4000，这与预期的答案非常接近。但这可能不是您的教授所期望的，因为他特别要求输出精确 -4000.4。

注意：显然，您可以使用另一种更精确的数字类型，但您的教授可能希望您使用double，因此我不会详细介绍。

编辑：正如评论中提到的@phuclv，一些编译器没有实现long double为 80 位浮点值，而是可能具有与double（64 位）相同的精度。因此，您可能需要寻找提供适当 80 位精度long double甚至128 位 IEEE-754 四倍精度浮点类型的库。虽然这绝对不会被认为是“容易”的。

（大部分是正确的）修复

好吧，你不能无限精确，因为double类型有 epsilon = 2^(-52)，但是你可以在加法中更聪明，而不仅仅是将大值添加到小值中（记住：大值“吸收”小值因为double浮点运算不精确）。基本上，您应该计算一个具有值的成对乘法的数组，然后对其进行排序（基于绝对值），然后使用std::accumulate以下方法添加值：

#include <iostream>
#include <numeric>
#include <vector>
#include <functional>
//Mind the use of these two new STL libraries
#include <algorithm> //std::sort and std::transform
#include <cmath> //abs()



int main(){

    std::vector<double> x{1.0e20, -1.0e3, 0.1, 1.0e20};
    std::vector<double> y{1.0, 4.0, -4.0, -1.0};
    //The vector with the pairwise products
    std::vector<double> products(x.size());

    //Do element-wise multiplication
    //C code: products[i] += x[i] * y[i];
    std::transform(x.begin(), x.end(), y.begin(), products.begin(), std::multiplies<double>());

    //Sort the array based on absolute-value
    auto sort_abs = [] (double a, double b) { return abs(a) < abs(b); };
    std::sort(products.begin(), products.end(), sort_abs);

    //Add the values of the products(note the init=0.0)
    double result = std::accumulate(products.begin(), products.end(), 0.0);

    std::cout << result << std::endl;
    return 0;
}

使用此新代码，结果如预期： -4000.4

坚韧它显然有它的局限性。例如，如果输入是向量 v1 = {100.0, 1E20} 和 v2 = {10.0, 1.0}，结果应该返回100000000000000001000，显然只会返回 1E20。