如何正确舍入平方根函数?
How to Correctly Round a Square Root Function?
我目前正在开发一个 Java 数学库,它将包括各种正确舍入的函数(即 sqrt、cbrt、exp、sin、gamma 和 ln)。我已经使用巴比伦方法编写了一个平方根算法,该算法与正确答案的误差在 1 ulp 以内。但是,我无法弄清楚如何正确计算应该以哪种方式舍入数字以表示对输入的实际平方根的最佳近似值。包含可以扩展到其他函数的原理的答案将是首选,但我听说 sqrt 比许多先验函数更简单,也将不胜感激专门的解决方案。
此外,这是我的代码在这个问题的原始提交时的清理版本:
public static double sqrt(double x) {
long bits = Double.doubleToLongBits(x);
// NaN and non-zero negatives:
if (Double.isNaN(x) || x < 0) return Double.NaN;
// +-0 and 1:
if (x == 0d || x == 1d) return x;
// Halving the exponent to come up with a good initial guess:
long exp = bits << 1;
exp = (exp - 0x7fe0000000000000L >> 1) + 0x7fe0000000000000L >>> 1 & 0x7ff0000000000000L;
double guess = Double.longBitsToDouble(bits & 0x800fffffffffffffL | exp);
double nextUp, nextDown, guessSq, nextUpSq, nextDownSq;
// Main loop:
while (true) {
guessSq = guess * guess;
if (guessSq == x) return guess;
nextUp = Math.nextUp(guess);
nextUpSq = nextUp * nextUp;
if (nextUpSq == x) return nextUp;
if (guessSq < x && x < nextUpSq) {
double z = x / nextUp;
if (z * nextUp > x) z = Math.nextDown(z);
return z < nextUp ? nextUp : guess;
}
nextDown = Math.nextDown(guess);
nextDownSq = nextDown * nextDown;
if (nextDownSq == x) return nextDown;
if (nextDownSq < x && x < guessSq) {
double z = x / guess;
if (z * guess > x) z = Math.nextDown(z);
return z < guess ? guess : nextDown;
}
// Babylonian method:
guess = 0.5 * (guess + x / guess);
}
}
如您所见,我正在使用除法作为测试。但是,我认为这需要除法向 0 舍入,这在 Java.
中显然不会发生。
根据泰勒定理,平方根函数局部近似为一个线性函数,斜率为1/2√x,为正。因此,您可以将误差与平方 x - (√x)² 中的误差相关联,其中 √x 被理解为近似根。然后你朝着最小化这个错误的方向四舍五入。
无论如何,x - (√x)² 的计算会发生灾难性的抵消,您可能需要更高的精度才能可靠地计算它。不确定这样做的好处是否值得。
我目前正在开发一个 Java 数学库,它将包括各种正确舍入的函数(即 sqrt、cbrt、exp、sin、gamma 和 ln)。我已经使用巴比伦方法编写了一个平方根算法,该算法与正确答案的误差在 1 ulp 以内。但是,我无法弄清楚如何正确计算应该以哪种方式舍入数字以表示对输入的实际平方根的最佳近似值。包含可以扩展到其他函数的原理的答案将是首选,但我听说 sqrt 比许多先验函数更简单,也将不胜感激专门的解决方案。
此外,这是我的代码在这个问题的原始提交时的清理版本:
public static double sqrt(double x) {
long bits = Double.doubleToLongBits(x);
// NaN and non-zero negatives:
if (Double.isNaN(x) || x < 0) return Double.NaN;
// +-0 and 1:
if (x == 0d || x == 1d) return x;
// Halving the exponent to come up with a good initial guess:
long exp = bits << 1;
exp = (exp - 0x7fe0000000000000L >> 1) + 0x7fe0000000000000L >>> 1 & 0x7ff0000000000000L;
double guess = Double.longBitsToDouble(bits & 0x800fffffffffffffL | exp);
double nextUp, nextDown, guessSq, nextUpSq, nextDownSq;
// Main loop:
while (true) {
guessSq = guess * guess;
if (guessSq == x) return guess;
nextUp = Math.nextUp(guess);
nextUpSq = nextUp * nextUp;
if (nextUpSq == x) return nextUp;
if (guessSq < x && x < nextUpSq) {
double z = x / nextUp;
if (z * nextUp > x) z = Math.nextDown(z);
return z < nextUp ? nextUp : guess;
}
nextDown = Math.nextDown(guess);
nextDownSq = nextDown * nextDown;
if (nextDownSq == x) return nextDown;
if (nextDownSq < x && x < guessSq) {
double z = x / guess;
if (z * guess > x) z = Math.nextDown(z);
return z < guess ? guess : nextDown;
}
// Babylonian method:
guess = 0.5 * (guess + x / guess);
}
}
如您所见,我正在使用除法作为测试。但是,我认为这需要除法向 0 舍入,这在 Java.
中显然不会发生。根据泰勒定理,平方根函数局部近似为一个线性函数,斜率为1/2√x,为正。因此,您可以将误差与平方 x - (√x)² 中的误差相关联,其中 √x 被理解为近似根。然后你朝着最小化这个错误的方向四舍五入。
无论如何,x - (√x)² 的计算会发生灾难性的抵消,您可能需要更高的精度才能可靠地计算它。不确定这样做的好处是否值得。