有效地计算（一 - K）/（A + K），具有改进的准确度 [英] Efficiently computing (a - K) / (a + K) with improved accuracy
问题描述
在各种情况下，例如，对于参数减少数学函数，一个需要计算（A - K）/（A + K），其中 A 是一个积极的可变参数和 K 是一个常数。在许多情况下， K 是二的幂，这是相关的我的工作中使用的情况下。我要寻找有效的方法来更准确地计算这个商数比使用简单的分工来完成。可以假设为乘加（FMA）硬件支持，因为这个操作是由所有主要的CPU和GPU架构在这个时候提供，并在C / C ++可通过功能 FMA（）和 fmaf（）。
In various contexts, for example for the argument reduction for mathematical functions, one needs to compute (a - K) / (a + K), where a is a positive variable argument and K is a constant. In many cases, K is a power of two, which is the use case relevant to my work. I am looking for efficient ways to compute this quotient more accurately than can be accomplished with the straightforward division. Hardware support for fused multiply-add (FMA) can be assumed, as this operation is provided by all major CPU and GPU architectures at this time, and is available in C/C++ via the functionsfma() and fmaf().
为了便于摸索，我与浮动算术试验。由于我打算端口的方法来双击算术为好，可以使用除了这两个参数和结果的本土precision使用更高的任何操作。我最好的解决办法，到目前为止是：
For ease of exploration, I am experimenting with float arithmetic. Since I plan to port the approach to double arithmetic as well, no operations using higher than the native precision of both argument and result may be used. My best solution so far is:
/* 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);
有关参数 A 在区间 [K / 2，4.23 * K] ，code以上计算商几乎是正确舍入的所有输入（最大误差是极为接近0.5 ULPS），前提是 K 是2的幂，并且没有溢或下溢中间结果。对于 K 不是一个两个电源，这code是仍比基于分工天真的算法更准确。在性能方面，这款code能的更快的比对平台上的幼稚的方法，把浮点倒数可以计算比浮点除法快。
For arguments a in the interval [K/2, 4.23*K], code above computes the quotient almost correctly rounded for all inputs (maximum error is exceedingly close to 0.5 ulps), provided that K is a power of 2, and there is no overflow or underflow in intermediate results. For K not a power of two, this code is still more accurate than the naive algorithm based on division. In terms of performance, this code can be faster than the naive approach on platforms where the floating-point reciprocal can be computed faster than the floating-point division.
我提出以下意见时， K = 2 N ：当上限的工作时间增加到 8 * K ， 16 * K ，...逐渐最大误差增加，并开始从下方慢慢接近幼稚计算的最大误差。不幸的是，同样不会出现为下界的间隔是真实的。如果下界滴 0.25 * K ，该改进的方法的最大误差以上等于幼稚方法的最大误差。
I make the following observation when K = 2n: When the upper bound of the work interval increases to 8*K, 16*K, ... maximum error increases gradually and starts to slowly approximate the maximum error of the naive computation from below. Unfortunately, the same does not appear to be true for the lower bound of the interval. If the lower bound drops to 0.25*K, the maximum error of the improved method above equals the maximum error of the naive method.
有没有计算Q A方法=（A - K）/（A + K），可以实现更小的最大错误（测 ULP VS的数学结果）相比，无论天真的方法和上面的code序列，在更宽的区间的特别是对区间的下限小于 0.5 * K ？的效率是重要，但多了一些操作要比在上述code使用可能可以容忍的。
Is there a method to compute q = (a - K) / (a + K) that can achieve smaller maximum error (measured in ulp vs the mathematical result) compared to both the naive method and the above code sequence, over a wider interval, in particular for intervals whose lower bound is less than 0.5*K? Efficiency is important, but a few more operations than are used in the above code can likely be tolerated.
在下面的一个答案，有人指出，我可以通过返回商作为两个操作数的未计算总和，也就是说，作为一个头尾对①提高精度：QLO ，即类似于著名的双 - 浮动和双 - 双击格式。在我的code以上，这将意味着改变最后一行 QLO = R * E 。
In one answer below, it was pointed out that I could enhance accuracy by returning the quotient as an unevaluated sum of two operands, that is, as a head-tail pair q:qlo, i.e. similar to the well-known double-float and double-double formats. In my code above, this would mean changing the last line to qlo = r * e.
这做法肯定是有用的，我已经想到了一个延伸期precision对数用于在使用 POW（）。但它并没有从根本上与在其上增强计算提供更准确的商数的时间间隔的期望加宽帮助。在我期待在一个特定的情况下，我想用 K = 2 （单precision）或 K = 4 （双precision），以保持主逼近区间狭窄，以及
This approach is certainly useful, and I had already contemplated its use for an extended-precision logarithm for use in pow(). But it doesn't fundamentally help with the desired widening of the interval on which the enhanced computation provides more accurate quotients. In a particular case I am looking at, I would like to use K=2 (for single precision) or K=4 (for double precision) to keep the primary approximation interval narrow, and the interval for a is roughly [0,28]. The practical problem I am facing is that for arguments < 0.25*K the accuracy of the improved division is not substantially better than with the naive method.
推荐答案
一种可能性是跟踪的m和p误差成m1和p1和古典德克尔/ Schewchuk：
One possibility is to track error of m and p into m1 and p1 with classical Dekker/Schewchuk:
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;
然后，更正天真的划分：
Then, correct the naive division:
q=m/p;
r0=fmaf(p,-q,m);
r1=fmaf(p1,-q,m1);
r=r0+r1;
q1=r/p;
q=q+q1;
这会花费你2师，而应该是接近一半ULP如果我没有搞砸了。
That'll cost you 2 divisions, but should be near half ulp if I didn't screw up.
但这些分歧可以用对逆乘法来代替没有任何问题，因为第一个错误四舍五入司将余数r进行补偿，二是不正确的四舍五入师并不真正的问题（修正Q1的最后一位赢得 T改变任何东西）。
But these divisions can be replaced by multiplications with inverse of p without any problem, since the first incorrectly rounded division will be compensated by remainder r, and second incorrectly rounded division does not really matter (the last bits of correction q1 won't change anything).
这篇关于有效地计算（一 - K）/（A + K），具有改进的准确度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
