algorithm - 如何修复此浮点平方根算法

Question

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

对于输入 0x7F7FFFFF (3.4028234663852886E38).，算法不会收敛到正确答案 18446743523953729536.000000 该算法的答案给出 18446743523953737728.000000。

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

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；这给了我们最坏的近似误差的另一个 half-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 并得到一个完全可表示的结果 --- s*s-foo 精确计算。

当您除以 s 时，您会引入额外的 half-ulp 误差，这里最多为 2^-48，因为我们的误差已经很小了。

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

要将其转化为正确舍入的保证，我们需要证明在 1/2 和 2 之间没有 float，除了我特例的那个，它的平方根是在两个连续 float 之间的中点的 3*2^-48 范围内。你可以做一些错误分析，得到一个丢番图方程，找到那个丢番图方程的所有解，找到它们对应的输入，并计算出算法对这些解的作用。 (如果你这样做，会有一个“物理”解决方案和一堆“非物理”解决方案。一个真正的解决方案是我唯一特例的解决方案。)然而，可能有更简洁的方法。