c++ - 执行双浮点除法的正确算法是什么？

Question

我正在遵循 this paper by Andrew Thall 提供的算法描述使用 df64 数据类型执行数学运算的算法，df64 数据类型是一对 32 位 float ，用于模拟 64 位 float 的精度。然而，他们编写除法和平方根函数的方式似乎存在一些不一致(错误？)。

论文中的除法函数是这样写的:

float2 df64_div(float2 B, float2 A) {
    float xn = 1.0f / A.x;
    float yn = B.x * xn;
    float diff = (df64_diff(B, df64_mult(A, yn))).x;
    float2 prod = twoProd(xn, diffTerm);

    return df64_add(yn, prodTerm);
}

^{用于编写此代码的语言似乎是 Cg，以供引用，但如果您将 float2 视为仅仅是一个别名，您应该能够用 C++ 解释此代码struct float2{float x, y;};，带有一些额外的语法来支持直接在类型上进行算术运算。}

作为引用，这些是此代码中使用的函数的 header :

float2 df64_add(float2 a, float2 b);
float2 df64_mult(float2 a, float2 b);
float2 df64_diff(/*Not provided...*/);
float2 twoProd(float a, float b);

因此，有几个问题立即凸显出来:

• diffTerm 和 prodTerm 从未定义。定义了两个变量，diff 和 prod，，但不确定这些是该代码中预期的术语。
• 未提供 df64_diff 声明。据推测，这是为了支持减法；但同样，这还不清楚。
• df64_mult 是一个不接受 32 位 float 作为参数的函数；它仅支持两对 32 位 float 作为参数。目前尚不清楚论文期望这个函数调用如何编译
• 对于 df64_add 也是如此，它也只接受成对的 32 位 float 作为参数，但在此处调用时第一个参数仅为单个 32 位 float 。

我有根据地猜测这是该代码的正确实现，但是因为即使该函数的正确实现在计算中也存在不可避免的错误，所以我无法判断它是否正确，即使它给出了值“似乎”正确:

float2 df64_div(float2 B, float2 A) {
    float xn = 1.0f / A.x;
    float yn = B.x * xn;
    float diff = (df64_diff(B, df64_mult(A, float2(yn, 0)))).x;
    float2 prod = twoProd(xn, diff);

    return df64_add(float2(yn, 0), prod);
}

float2 df64_diff(float2 a, float2 b) {
    return df64_add(a, float2(-b.x, -b.y));
}

所以我的问题是:论文中看到的这个算法的书面实现是否准确(因为它取决于我不知道的 Cg 语言的行为？)，或者不是？不管怎样，我对该代码的插值是否是论文中描述的除法算法的正确实现？

注意:我的目标语言是 C++，因此虽然语言之间的差异(对于这种算法)很小，但我的代码是用 C++ 编写的，并且我正在寻找 C++ 语言的正确性。

Answer 1

回顾书中所写的伪代码算法似乎支持该算法的 C++ 实现，尽管我对 Cg 不熟悉意味着我无法证明该实现对于 Cg 是正确的。

因此，将这些步骤分解为简单的英语:

1. 该函数采用两个参数，每个参数都是[伪] double 浮点值，并且第二个参数不等于 0
2. 变量 x_n 被赋予[伪]双除数的高阶分量的算术倒数，使用​​单精度 float 学计算
3. 变量 y_n 被赋予[伪]双倍被除数的高阶分量与 x_n 的乘积，使用单精度 float 学计算
4. 计算[伪]双除数和 y_n 的乘积
   • 这是第一个棘手的部分，因为本文没有描述[伪]双 x 单乘法的算法。我们可以在Cg算法中看到，Cg算法清楚地映射到这一步1对1，但是将标量值提升为 vector 值的Cg规则是未知的。
   • 但是，我们可以说的是，我们确实有一个将 double 乘以 double 的函数，并且可以通过用 0 填充其低阶分量来将单精度数提升为 double ，因此我们可以做到这一点。
5. 计算 Dividend 与步骤 4 中计算出的乘积之间的差值，仅将高阶分量保留为单精度浮点值
   • 让这个问题变得棘手的是，这篇论文没有描述减法算法。然而，它确实描述了一种将 [IEEE754-]double 转换为 [pseudo-]double 的算法，我们可以观察到负的 [IEEE754-]double 在转换时，其高阶和高阶都具有负值。低阶分量。因此从逻辑上讲，[伪]double 可以通过简单地否定它的两个组成部分来否定。添加一个负数在数学上相当于减法，因此我们可以利用这些知识构建减法算法。
6. 执行 x_n 和步骤 5 的乘积，保留扩展精度，否则会在单个 x 乘法中丢失。
   • twoProd函数的存在正是为了这个目的。
7. 计算第 6 步与 y_n 的总和
   • 同样，如果我们简单地通过用 0 填充低阶分量将 y_n 提升为 [pseudo-]double，我们就可以使用 [pseudo-]double 加法算法
8. 第7步的结果就是返回值

因此，了解这个算法后，我们可以将每个步骤直接映射到我编写的 C++ 算法:

//(1) Takes two [pseudo-]doubles, returns a [pseudo-]double
float2 df64_div(float2 B, float2 A) {
    //(2) single float divided by single float
    float xn = 1.0f / A.x;
    // (3) single float multiplied by single float
    float yn = B.x * xn;
    //                        (4) double x double multiplication
    //                                       (4a) yn promoted to [pseudo-]double
    //            (5) subtraction                           (5a) only higher order component kept
    float diff = (df64_diff(B, df64_mult(A, float2(yn, 0)))).x;
    // (6) single x single multiplication with extra precision preserved using twoProd
    float2 prod = twoProd(xn, diff);
    // (7) adding higher-order division to lower order division
    //              (7a) yn promoted to [pseudo-]double
    // (8) value is returned
    return df64_add(float2(yn, 0), prod);
}

float2 df64_diff(float2 a, float2 b) {
    //                 (5a) negating both components is a logical negation of the whole number
    return df64_add(a, float2(-b.x, -b.y));
}

由此，我们可以得出结论，这是本文中描述的算法的正确实现，我所做的一些测试证实了以这种方式执行这些操作会产生看似正确的结果。