如何修复这个浮点平方根算法

Question

我正在尝试计算各种输入的 IEEE-754 32 位浮点平方根，但对于一个特定输入，以下基于 Newton-Raphson 方法的算法不会收敛，我想知道我能做什么解决问题？对于我正在设计的平台，我有一个 32 位浮点加法器/减法器、乘法器和除法器。

对于输入 0x7F7FFFFF (3.4028234663852886E38)。，算法不会收敛到正确答案 18446743523953729536.000000 这个算法的答案是 18446743523953737728.000000。

在硬件实现之前，我使用 MATLAB 来实现我的代码。我只能使用单精度浮点值，（所以没有双精度）。

clc; clear; close all;

% Input
R = typecast(uint32(hex2dec(num2str(dec2hex(((hex2dec('7F7FFFFF'))))))),'single')

% Initial estimate
OneOverRoot2 = single(1/sqrt(2));
Root2 = single(sqrt(2));

% Get low and high bits of input R
hexdata_high = bitand(bitshift(hex2dec(num2hex(single(R))),-16),hex2dec('ffff'));
hexdata_low = bitand(hex2dec(num2hex(single(R))),hex2dec('ffff'));

% Change exponent of input to -1 to get Mantissa
temp = bitand(hexdata_high,hex2dec('807F'));
Expo = bitshift(bitand(hexdata_high,hex2dec('7F80')),-7);
hexdata_high = bitor(temp,hex2dec('3F00'));
b = typecast(uint32(hex2dec(num2str(dec2hex(((bitshift(hexdata_high,16)+ hexdata_low)))))),'single');

% If exponent is odd ...
if (bitand(Expo,1))
    % Pretend the mantissa [0.5 ... 1.0) is multiplied by 2 as Expo is odd,
    %   so it now has the value [1.0 ... 2.0)
    % Estimate the sqrt(mantissa) as [1.0 ... sqrt(2))
    % IOW: linearly map (0.5 ... 1.0) to (1.0 ... sqrt(2))
    Mantissa = (Root2 - 1.0)/(1.0 - 0.5)*(b - 0.5) + 1.0;
else
    % The mantissa is in range [0.5 ... 1.0)
    % Estimate the sqrt(mantissa) as [1/sqrt(2) ... 1.0)
    % IOW: linearly map (0.5 ... 1.0) to (1/sqrt(2) ... 1.0)
    Mantissa = (1.0 - OneOverRoot2)/(1.0 - 0.5)*(b - 0.5) + OneOverRoot2;
end

newS = Mantissa*2^(bitshift(Expo-127,-1));
S=newS

% S = (S + R/S)/2 method
for j = 1:6 
    fprintf('S  %u %f %f\n', j, S, (S-sqrt(R)));
    S = single((single(S) + single(single(R)/single(S))))/2;
    S = single(S);
end

goodaccuracy =  (abs((single(S)-single(sqrt(single(R)))))) < 2^-23
difference = (abs((single(S)-single(sqrt(single(R))))))

% Get hexadecimal output
hexdata_high = (bitand(bitshift(hex2dec(num2hex(single(S))),-16),hex2dec('ffff')));
hexdata_low = (bitand(hex2dec(num2hex(single(S))),hex2dec('ffff')));
fprintf('FLOAT: T  Input: %e\t\tCorrect: %e\t\tMy answer: %e\n', R, sqrt(R), S);
fprintf('output hex = 0x%04X%04X\n',hexdata_high,hexdata_low);
out = hex2dec(num2hex(single(S)));

Answer 1

我对此感到震惊。这是我想出的：

float mysqrtf(float f) {
  if (f < 0) return 0.0f/0.0f;
  if (f == 1.0f / 0.0f) return f;
  if (f != f) return f;

  // half-ass an initial guess of 1.0.
  int expo;
  float foo = frexpf(f, &expo);
  float s = 1.0;
  if (expo & 1) foo *= 2, expo--;

  // this is the only case for which what's below fails.
  if (foo == 0x0.ffffffp+0) return ldexpf(0x0.ffffffp+0, expo/2);

  // do four newton iterations.
  for (int i = 0; i < 4; i++) {
   float diff = s*s-foo;
    diff /= s;
    s -= diff/2;
  }

  // do one last newton iteration, computing s*s-foo exactly.
  float scal = s >= 1 ? 4096 : 2048;
  float shi = (s + scal) - scal; // high 12 bits of significand
  float slo = s - shi; // rest of significand
  float diff = shi * shi - foo; // subtraction exact by sterbenz's theorem
  diff += 2 * shi * slo; // opposite signs; exact by sterbenz's theorem
  diff += slo * slo;
  diff /= s; // diff == fma(s, s, -foo) / s.
  s -= diff/2;

  return ldexpf(s, expo/2);
}

首先要分析的是浮点运算中的公式(s*s-foo)/s。如果s 是sqrt(foo) 的一个足够好的近似值，Sterbenz 定理告诉我们分子在正确答案的 ulp(foo) 范围内 --- 所有这些误差都是计算 s*s 的近似误差。然后我们除以s；这在最坏的情况下给了我们另一个半 ulp 的近似误差。因此，即使没有融合乘加，diff 也在应有的 1.5 ulp 范围内。然后我们将其除以二。

请注意，最初的猜测本身并不重要，只要你用足够的牛顿迭代来跟进它。

用 abs(s - foo/s) 测量近似值 s 到 sqrt(foo) 的误差。我最初猜测的 1 的误差最多为 1。精确算术中的牛顿迭代将误差平方并除以 4。浮点算术中的牛顿迭代——我做四次的那种——平方错误，将其除以 4，然后再产生 0.75 ulp 的错误。你这样做了四次，你发现你最多有一个相对误差0x0.000000C4018384，大约是 0.77 ulp。这意味着四次牛顿迭代产生了一个忠实的四舍五入结果。

我做了第五个牛顿步骤来得到一个正确舍入的平方根。它起作用的原因有点复杂。

shi 持有s 的“上半部分”，而slo 持有“下半部分”。每个有效数字的最后 12 位将为零。这尤其意味着shi * shi 和shi * slo 和slo * slo 完全可以表示为floats。

s*s 与 foo 相差两个 ulps。 shi*shi 在 s*s 的 2047 ulps 内。因此shi * shi - foo 在零的 2049 ulps 内；特别是，它是完全可表示的并且小于 2^-10。

您可以检查您是否可以添加 2 * shi * slo 并获得在零的 2^-22 内的完全可表示的结果，然后添加 slo*slo 并获得完全可表示的结果 --- @ 987654343@ 计算准确。

当你除以s 时，你会额外增加半个 ulp 的错误，这里最多为 2^-48，因为我们的错误已经很小了。

现在我们进行牛顿步。我们已经正确计算出当前误差在 2^-46 以内。将其中的一半添加到 s 中，我们得到的平方根在 3*2^-48 以内。

为了保证正确舍入，我们需要证明在1/2和2之间不存在floats，除了我特例的那个，它的平方根在3*2以内-48 两个连续的floats 之间的中点。你可以做一些错误分析，得到一个丢番图方程，找到这个丢番图方程的所有解，找到它们对应的输入，然后计算出算法对这些的作用。 （如果你这样做，就会有一个“物理”解决方案和一堆“非物理”解决方案。一个真正的解决方案是我唯一特殊的解决方案。）但是，可能有一种更清洁的方法。