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

我试图计算各种输入的IEEE-754 32位浮点平方根，但是对于一个特定的输入，下面基于Newton-Raphson方法的算法不会收敛，我想知道我能做些什么来解决这个问题对于我正在设计的平台，我有一个32位浮点加法器/减法器、乘法器和除法器。
对于输入0x7f7fffff（3.4028234663852886e38），算法不会收敛到1844674352395375336.000000的正确答案，该算法的答案给出1844674352395737728.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)));

最佳答案

我猛击了一下我想到的是：

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);
}

首先要分析的是浮点运算中的公式如果CC是一个足够好的近似到cc，斯特本兹的定理告诉我们分子在正确答案的一个（）中，所有这些误差都是计算的误差。然后我们用“CC”来划分，这给了我们最差的另一半误差。因此，即使没有熔合乘法加法，(s*s-foo)/s也在1.5ulp之内。我们把它除以2。
注意，只要你用足够多的牛顿迭代来跟踪它，最初的猜测本身并不重要。
用ABS（S—Fo/S）测量近似S到SqRT（FoO）的误差。我最初猜测1的误差最多是1。精确算术中的牛顿迭代将误差平方并除以4在浮点运算中，牛顿迭代法——我做了四次——将误差平方，除以4，再产生0.75ULP的误差你这样做了四次，你发现你最多有一个相对误差s，大约是0.77ulp这意味着四次牛顿迭代得到一个精确的四舍五入结果。
我做了第五步牛顿得到一个正确的四舍五入平方根。它工作的原因有点复杂。
sqrt(foo)保留s*s的“上半部分”，而s保留“下半部分”。每个有效位中的最后12位将为零。这特别意味着diff和0x0.000000C4018384以及shi可以精确地表示为ss。
slo在shi * shi的两个ulps内shi * slo在slo * slo的2047 ulps内因此float在2049ulps以内；特别是，它完全可以表示并且小于2-10。
您可以检查是否可以添加s*s并获得精确可表示的结果（在0的2-22范围内），然后添加foo并获得精确可表示的结果---shi*shi计算结果。
当你除以s*s时，你又会产生一个额外的半ULP误差，因为我们的误差已经很小了，所以这里最多是2-48。
现在我们做牛顿步。我们已经计算出当前的误差在2-46之间。把它的一半加到shi * shi - foo中，我们得到的平方根在3*2-48范围内。
为了保证舍入正确，我们需要证明在1/2到2之间没有2 * shi * slos，除了一个特殊的cased，它的平方根在两个连续的slo*slos之间中点的3*2-48之内。你可以做一些误差分析，得到一个丢番图方程，找到丢番图方程的所有解找出它们对应的输入，并计算出算法在这些输入上的作用。（如果你这样做，有一个“物理”解决方案和一堆“非物理”解决方案唯一真正的解决方案是我唯一的特例。）不过，可能有一个更干净的方法。