恐怕我不知道如何保证正确舍入的双精度多维数据集根，但是可以提供一个在问题中也提到过几乎正确舍入的根。换句话说，最大误差似乎非常接近0.5 ulp。

彼得·马克斯坦（Peter Markstein），“ IA-64和基本功能：速度和精度”（Prentice-Hall，2000年）

提出了基于FMA的有效技术，可以正确舍入倒数，平方根和倒数平方根，但是在这方面它不涵盖立方根。通常，Markstein的方法需要一个初步的结果，该结果要精确到最终舍入序列之前的1
ulp以内。我没有数学上的手段将他的技术扩展到立方根的修整，但是在我看来，从原理上讲这应该是可能的，这是与倒数平方根相似的挑战。

按位算法很容易通过正确的舍入来计算根。由于不会发生IEEE-754舍入为最接近或什至平均的舍入情况，因此只需要进行计算，直到产生所有尾数位加一个舍入位即可。基于二项式定理的平方根逐位算法在非还原和还原变体中都是众所周知的，并且已成为硬件实现的基础。通过二项式定理的相同方法适用于立方根，并且鲜为人知的论文提出了非还原实现的详细信息：

H. Peng，“提取平方根和立方根的算法”，第5届IEEE国际计算机算术研讨会论文集，第121-126页，1981年。

从我的实验中可以看出，这对于从整数提取立方根非常有效。由于每次迭代仅产生一个结果位，因此速度不是很快。对于浮点运算中的应用，其缺点是使用了两个簿记变量，它们需要最终结果位数的大约两倍。这意味着需要使用128位整数算术来实现双精度立方根。

我下面的C99代码基于哈雷针对立方根的有理方法，该方法具有三次收敛性，这意味着初始逼近不必非常精确，因为有效位数在每次迭代中都增加了三倍。可以以各种方式来安排计算。通常，将迭代方案安排为

new_guess：= old_guess +更正

因为对于足够接近的初始猜测，correction它明显小于old_guess。这导致多维数据集根的以下迭代方案：

x：= x-x *（x 3 -a）/（2 * x 3 + a）

这种特殊的安排也列在Kahan关于cube
root的注释中。它具有进一步自然地将自身借给FMA（融合乘加）的优点。一个缺点是2 * x
3的计算可能会导致溢出，因此对于双精度输入域的至少一部分需要使用自变量约简方案。在我的代码中，我基于对IEEE-754双精度操作数的指数的直接操作，将参数归约简单地应用于
所有 非异常输入。

间隔[0.125,1）作为一次近似间隔。使用多项式最小极大逼近，可返回[0.5,1]中的初始猜测。较窄的范围便于将单精度算法用于计算的低精度部分。

我无法证明实现的错误范围，但是，使用2亿个随机测试向量针对参考实现（精确到200位左右）进行测试，发现总共有277个错误的舍入结果（错误率约为1.4ppm）最大误差为0.500012 ulps。

double my_cbrt (double a)
{
    double b, u, v, r;
    float bb, uu, vv;
    int e, f, s;

    if ((a == 0.0) || isinf(a) || isnan(a)) {
        /* handle special cases */
        r = a + a;
    } else {
        /* strip off sign-bit */
        b = fabs (a);
        /* compute exponent adjustments */
        b = frexp (b, &e);
        s = e - 3*342;
        f = s / 3;
        s = s - 3 * f;
        f = f + 342;
        /* map argument into the primary approximation interval [0.125,1) */
        b = ldexp (b, s);
        bb = (float)b;
        /* approximate cube root in [0.125,1) with relative error 5.22e-3 */
        uu =                0x1.2f32c0p-1f;
        uu = fmaf (uu, bb, -0x1.62cc2ap+0f);
        uu = fmaf (uu, bb,  0x1.7546e0p+0f);
        uu = fmaf (uu, bb,  0x1.5d0590p-2f);
        /* refine cube root using two Halley iterations w/ cubic convergence */
        vv = uu * uu;
        uu = fmaf (fmaf (vv, uu, -bb) / fmaf (vv, 2.0f*uu, bb), -uu, uu);
        u = (double)uu;
        v = u * u; // this product is exact
        r = fma (fma (v, u, -b) / fma (v, 2.0*u, b), -u, u);
        /* map back from primary approximation interval by jamming exponent */
        r = ldexp (r, f);
        /* restore sign bit */
        r = copysign (r, a);
    }
    return r;
}