logo标签列表

高效计算 (a - K) / (a + K) 并提高精度

calgorithmfloating-pointfloating-accuracy
21
在不同的场景中，例如对于数学函数参数的减少，需要计算 (a - K) / (a + K)，其中 a 是一个正的变量参数，而 K 是一个常数。 在许多情况下，K 是二的幂，这是我的工作相关用例。我正在寻找比直接除法更有效的方法来更准确地计算这个商。可以假设硬件支持融合乘加（FMA），因为此操作由所有主要的 CPU 和 GPU 架构提供，在 C/C++ 中可通过函数 fma() 和 fmaf()获得。
为了便于探索，我正在尝试使用 float 算法。由于我计划将方法移植到 double 算法中，因此不能使用高于参数和结果本机精度的操作。目前为止，我最好的解决方案是：
 /* Compute q = (a - K) / (a + K) with improved accuracy. Variant 1 */
 m = a - K;
 p = a + K;
 r = 1.0f / p;
 q = m * r;
 t = fmaf (q, -2.0f*K, m);
 e = fmaf (q, -m, t);
 q = fmaf (r, e, q);

对于区间[K/2, 4.23*K]内的参数a，上述代码计算出所有输入的商几乎正确舍入（最大误差极接近0.5 ulps），前提是K为2的幂，并且中间结果没有溢出或下溢。对于非2的幂的K，此代码仍比基于除法的朴素算法更精确。在性能方面，在浮点倒数可以比浮点除法更快地计算的平台上，此代码可能比朴素方法更快。
当K=2n时，我做出以下观察：当工作区间的上限增加到8*K、16*K等时，最大误差逐渐增加，并开始从下方缓慢逼近朴素计算的最大误差。不幸的是，对于区间的下限却不是这样。如果下限降至0.25*K，改进后的方法的最大误差等于朴素方法的最大误差。
是否有一种方法来计算q=(a-K)/(a+K)，以实现与朴素方法和上述代码序列相比，更小的最大误差（用ulp来衡量）并且适用于更广的区间，特别是对于下限小于0.5*K的区间？效率很重要，但比上述代码中使用的操作稍多一些可能是可以容忍的。
在下面的一个答案中，指出我可以通过将商作为两个操作数的未求值和返回，即作为头尾对q:qlo来提高精度，类似于众所周知的双float和双double格式。在我的上述代码中，这意味着将最后一行改为qlo=r*e。
这种方法确实有用，我已经考虑过在pow()中使用它来扩展精度对数。但它并不能从根本上帮助扩大增强计算提供更准确商的区间。在我正在研究的一个特定情况中，我想使用K=2（单精度）或K=4（双精度）来保持主要近似区间较窄，而a的区间大约为[0,28]。我面临的实际问题是，对于小于0.25*K的参数，改进后的除法的精度与朴素方法的精度没有太大差别。
我不确定你所说的“平均误差曲线”是什么意思。我有兴趣将以ulps为单位测量的最大误差最小化。我通过对测试区间进行全面测试来确定误差，这就是为什么我在探索性工作中使用单精度算术的原因。 - njuffa
2
我想知道是否值得关注以下相对误差：(a / (a + k)) - (k / (a + k))？ - Brett Hale
@BrettHale 以这种方式重写表达式会导致最大ulp误差爆炸，因为当a接近K时会发生减法抵消。 - njuffa
如果双精度不是过于昂贵的话，那么在高精度下运行快速近似并提前终止是有意义的。当K是2的幂时，这非常可行（因为除以K很便宜）。 - user3528438
1
不幸的是，在某些平台上，“double”操作要昂贵得多（高达“float”操作的32倍）。由于我还想在“double”上使用相同的算法，因此没有便宜的“四倍”操作可用。因此，只需使用“本地”宽度操作（这也使向量化更容易）即可满足要求。 - njuffa
显示剩余4条评论
6个回答
4
如果a相对于K很大，那么(a-K)/(a+K) = 1 - 2K / (a + K)将会给出一个很好的近似值。如果a相对于K很小，那么2a / (a + K) - 1将会给出一个很好的近似值。如果K/2 ≤ a ≤ 2K，那么a-K是一个精确的操作，所以进行除法运算将会得到一个不错的结果。
如果您能够建议三条建议代码路径之间的切换点，我将很高兴通过我的测试框架运行它。虽然多分支代码不一定友好于向量化，因此可能效率低下，但在这种情况下，该问题可以通过预测来解决。 - njuffa
抱歉，我忽略了切换点已经足够指定的事实。我将算法翻译成了如下的C代码，并发现在[0.5K,4K)上的最大ulp误差仅略低于2.5 ulps，这比朴素方法更大：m = a - K; p = a + K; if ((0.5f*K <= a) && (a <= 2.0f*K)) { q = m / p; } else if (a < 0.5f*K) { q = 1.0f - 2.0f*K / p; } else { q = (2.0f * a) / p - 1.0f; } - njuffa
4

一种可能的方法是使用经典的Dekker/Schewchuk方法将m和p的误差跟踪到m1和p1中：

m=a-k;
k0=a-m;
a0=k0+m;
k1=k0-k;
a1=a-a0;
m1=a1+k1;

p=a+k;
k0=p-a;
a0=p-k0;
k1=k-k0;
a1=a-a0;
p1=a1+k1;

接下来，纠正这种幼稚的分割：

q=m/p;
r0=fmaf(p,-q,m);
r1=fmaf(p1,-q,m1);
r=r0+r1;
q1=r/p;
q=q+q1;

这将花费你2个除以2的部分，但如果我没搞砸的话，应该接近半个ulp。

但是，这些除法可以用p的倒数进行乘法替代，而不会有任何问题，因为第一个错误舍入的除法将被余数r补偿，第二个错误舍入的除法实际上并不重要（校正q1的最后几位不会改变任何东西）。

1
这似乎基本上是Simon Byrne建议的div2方法（https://dev59.com/J1sW5IYBdhLWcg3wFT2n#35434301），包括两个除法在内的18个操作。这已经完全编码了。我的实验表明，在[0.5 * K，32 * K）上，最大误差非常接近0.5 ulp，因此当增加区间的上限时，这似乎做得很好。然而，将下限降低到0.25 * K会将最大ulp误差增加到略小于2 ulps，比naive方法的最大误差约为1.625 ulp更糟糕。这可修复吗？ - njuffa
啊，看起来我把错误m1的符号搞错了...让我再检查一下。现在我编辑了我的答案，应该会更好。 - aka.nice
借助FMA的帮助，可以编写一个双精度浮点数除法，只需要进行一次倒数运算，而不是两次完整的除法。我怀疑在这里也可能有类似的优化。 - njuffa
3

我并没有一个确切的答案（精确的浮点数误差分析非常繁琐），但有几点观察：

  • 快速倒数指令（例如 RCPSS）不如除法准确，因此如果使用这些指令，您可能会看到准确度降低。
  • 如果 a ∈ [0.5×Kb, 21+n×Kb)，其中 Kb 是 K 下面的 2 的幂次方（如果 K 是 2 的幂次方，则为 K 本身），n 是 K 的有效数字中尾随零的数量（即如果 K 是 2 的幂次方，则 n=23），则可以精确计算 m。
  • 这类似于 Dekker（1971） 的简化版 div2 算法：要扩展范围（特别是下限），您可能需要从中加入更多的校正项（即将 m 存储为 2 个 float 的和，或使用一个 double）。
2
我熟悉关于快速求倒数的权衡。通常，硬件指令和适当数量的 NR 步骤的组合可以获得一个几乎完全舍入的倒数，即最大误差非常接近 0.5 ulps，这是可行的。在其他平台上，使用适当的除法加上一些 FMA 的相对较小开销在性能方面仍然是可以接受的。我知道 Dekker 的工作，但基本上只使用了其中的加法和乘法部分。我会再看一下，看看 div2 是否适用。 - njuffa
1
你是正确的：快速倒数不会有太大的差别，因为有修正项。 - Simon Byrne
我研究了一下双精度float除法，看起来至少需要13个操作。如果我只需要一个float结果，我可以节省两个操作。但是，我需要至少6个额外的操作来计算a+K和a-K，因此这种方法需要至少17个操作，而我的当前代码只需要7个操作。这似乎是最后的备选方案，性能影响很难证明。 - njuffa
我根据使用双float算术计算所有中间计算的方法编写了代码。不幸的是，我需要11个操作来计算作为两个双float运算对象的a+K和a-K。然后对它们进行除法需要11个操作，仅需要一个倒数，总共需要22个操作，比问题中使用7个操作的代码多15个操作。对于快速测试，我选择了区间[K/128，128*K)，这非常好，最大误差非常接近0.5 ulp。 - njuffa
2

因为我的目标只是扩大准确结果的间隔，而不是找到适用于所有可能值的解决方案，因此，在所有中间计算中使用双倍-float运算过于昂贵。

再思考一下这个问题，很明显，除法余数的计算，即我的问题中的e，是实现更精确结果的关键部分。在数学上，余数是(a-K) - q * (a+K)。在我的代码中，我只是使用m来表示(a-K)，并将(a+k)表示为m + 2*K，因为这样可以得到数值上优越的结果。

通过相对较小的额外计算成本，(a+K)可以表示为双倍-float，即头尾对p:plo，这导致了我的原始代码的以下修改：

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 2 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
mx = fmaxf (a, K);
mn = fminf (a, K);
plo = (mx - p) + mn;
t = fmaf (q, -p, m);
e = fmaf (q, -plo, t);
q = fmaf (r, e, q);

测试表明，在[a/2，224*a)范围内，这几乎可以提供正确舍入的结果，从而允许将准确结果实现的上限大幅增加。

将下限扩大到更精确地表示(a-K)需要。我们可以将其计算为双精度浮点头尾对m:mlo，这导致以下代码变体:

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 3 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
plo = (a < K) ? ((K - p) + a) : ((a - p) + K);
mlo = (a < K) ? (a - (K + m)) : ((a - m) - K);
t = fmaf (q, -p, m);
e = fmaf (q, -plo, t);
e = e + mlo;
q = fmaf (r, e, q);

详尽的测试表明，这个方法在区间[K/224, K*224)内能够提供几乎正确舍入的结果，但与我问题中的代码相比，需要多付出十次额外的操作，这是一个很高的代价，只是为了将最大误差从约1.625 ulps的天真计算降至接近0.5 ulp。

与我在问题中的原始代码一样，我们可以用(a-K)来表示(a+K)，从而消除对p和plo的尾部计算。这种方法得到的代码如下：

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 4 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
mlo = (a < K) ? (a - (K + m)) : ((a - m) - K);
t = fmaf (q, -2.0f*K, m);
t = fmaf (q, -m, t);
e = fmaf (q - 1.0f, -mlo, t);
q = fmaf (r, e, q);

如果主要关注的是将区间下限降低，这对于我提出的问题非常有帮助。对单精度情况进行全面测试表明，当K=2n时，在区间[a code=" "]K/224[/a], [a code=" "]4.23*K[/a]中的值都能够产生几乎正确的结果。这需要进行14或15次操作（取决于架构是否支持完整预测或者只是条件移动），比我的原始代码多出七到八个操作。

最后，为了避免计算m和p时固有的误差，可以直接将剩余计算基于原始变量a。对于K = 2n，以下代码可以计算在区间[a code=" "]K/224[/a], [a code=" "]K/3[/a])中的a值，产生几乎正确的结果：

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 5 */
m = a - K;
p = a + K;
r = 1.0f / p;       
q = m * r;
t = fmaf (q + 1.0f, -K, a);
e = fmaf (q, -a, t);
q = fmaf (r, e, q);
1
问题在于 (a + K) 的加法运算。任何在 (a + K) 中的精度损失都会被除法放大。问题并不在于除法本身。
如果 a 和 K 的指数相同（几乎相同），则不会丢失精度，如果指数之间的绝对差大于有效数字大小，则 (a + K) == a（如果 a 的数量级更大）或 (a + K) == K（如果 K 的数量级更大）。
没有办法阻止这种情况。增加有效数字位数（例如，在80x86上使用80位“扩展双精度”）只会稍微扩大“准确结果范围”。要理解为什么，请考虑smallest + largest（其中smallest是32位浮点数可以是最小正规的值）。在这种情况下（对于32位浮点数），您需要大约260位的有效数字位数才能完全避免精度损失。即使执行temp = 1/(a + K); result = a * temp - K / temp;也无济于事，因为您仍然面临完全相同的(a + K)问题（但它可以避免(a - K)中的类似问题）。此外，您不能执行result = anything / p + anything_error/p_error，因为除法不起作用。

我只能想到3种选择来接近所有可能适合32位浮点数的正值的0.5 ulps。但都不太可行。

第一种选择是预先计算每个值的查找表（使用“大实数”数学），对于32位浮点数，这将占用约2 GiB的空间（对于64位浮点数则完全疯狂）。当然，如果可能的a值范围小于“任何可以适合32位浮点数的正值”，则查找表的大小将被减小。
第二种选择是在运行时使用其他东西（“大实数”）进行计算，并转换为/从32位浮点数。
第三种选择涉及“某些东西”（我不知道它叫什么，但它很昂贵）。将舍入模式设置为“向正无穷舍入”，并计算temp1 = (a + K); if(a < K) temp2 = (a - K);然后切换到“向负无穷舍入”并计算if(a >= K) temp2 = (a - K); lower_bound = temp2 / temp1;。接下来执行a_lower = a并尽可能减少a_lower的值，并重复“lower_bound”计算，直到获得lower_bound的不同值，然后恢复到先前的a_lower值。之后，您基本上要执行相同的操作（但舍入模式相反，并且是增加而不是减少），以确定upper_bound和a_upper（从原始值a开始）。最后，像这样插值：a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range;。请注意，您将希望计算出初始的上限和下限，并在它们相等时跳过所有内容。还要警告你，这完全是“理论上，没有经过测试”的，我可能在某个地方搞砸了。

我想说的主要是（在我的观点下），你应该放弃并接受你无法达到0.5 ULP。抱歉.. :)

1
如果您可以放宽 API 返回另一个变量以模拟错误，那么解决方案就会变得简单得多：
float foo(float a, float k, float *res)
{
    float ret=(a-k)/(a+k);
    *res = fmaf(-ret,a+k,a-k)/(a+k);
    return ret;
}

这个解决方案只处理除法的截断误差，但不处理 a+k 和 a-k 的精度损失。

为了处理这些错误，我认为需要使用双精度或 bithack 来使用定点数。

测试代码已更新以人工生成输入中的非零最低有效位。

测试代码

https://ideone.com/bHxAg8

我猜你所说的“用其他变量来模拟误差”基本上是指将商作为头尾对（双浮点数，双倍精度）返回？我可以很容易地做到这一点（在我的代码中，这意味着用qlo = r * e替换最后一行），但我不知道它如何解决当下限边界低于0.5*K时误差急剧增加的问题。除法通常在任何平台上都很昂贵，我希望避免进行两次除法；倒数后跟随两个背向乘积可以获得更好的性能，因此我使用了这种方法。我会查看你的代码以探索细节。 - njuffa
我的测试框架通过对区间[0.5K, 4K)进行详尽的测试表明，上述代码计算商（视为未求和的ret:res）的最大误差不到1 ulp，这比朴素计算（约1.62 ulps）要好，但不如我问题中的代码（接近0.5 ulp）。我使用K=2进行测试，但只要没有下溢/上溢，任何2的幂都应该同样有效。如果您的测试结果与我的有实质性差异，请告诉我。 - njuffa
@njuffa 不，我同意您的测试结果。这也是我早些时候删除了这个答案的原因，因为我认为它并不能很好地解决问题。 - user3528438
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接