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

Question

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

Division函数在论文中是这样写的：

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 视为只是 struct float2{float x, y;}; 的别名，您应该能够在 C++ 中解释此代码一些额外的语法来支持直接在类型上的算术运算。}

作为参考，这些是此代码中使用的函数的标头：

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 位浮点数作为参数的函数；它只支持两对 32 位浮点数作为参数。目前尚不清楚论文期望这个函数调用如何编译
• df64_add 也一样，它也只接受成对的 32 位浮点数作为参数，但在这里调用的第一个参数只有一个 32 位浮点数。

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

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

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

p>

Answer 2

Xirema's answer 将 Thall 的高基长手除法算法忠实地呈现为 C++。基于对高精度参考的相当广泛的测试，我发现它的最大相对误差约为 2^-45，前提是中间计算中没有下溢。

在提供融合乘加运算 (FMA) 的平台上，由于Nagai et. al.，以下基于 Newton-Raphson 的除法算法可能更有效，并且在我的测试中达到相同的精度，即最大相对误差为 2^-45.

/*
  T. Nagai, H. Yoshida, H. Kuroda, Y. Kanada, "Fast Quadruple Precision 
  Arithmetic Library on Parallel Computer SR11000/J2." In: Proceedings 
  of the 8th International Conference on Computational Science, ICCS '08, 
  Part I, pp. 446-455.
*/
float2 div_df64 (float2 a, float2 b)
{
    float2 t, c;
    float r, s;
    r = 1.0f / b.x;
    t.x = a.x * r;
    s = fmaf (-b.x, t.x, a.x);
    t.x = fmaf (r, s, t.x);
    t.y = fmaf (-b.x, t.x, a.x);
    t.y = a.y + t.y;
    t.y = fmaf (-b.y, t.x, t.y);
    s = r * t.y;
    t.y = fmaf (-b.x, s, t.y);
    t.y = fmaf (r, t.y, s);
    c.x = t.x + t.y;
    c.y = (t.x - c.x) + t.y;
    return c;
}