2

我正在尝试计算各种输入的 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)));
4

1 回答 1

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 * shiandshi * slo和slo * slo可以精确地表示为floats。

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

您可以检查是否可以添加2 * shi * slo并获得在零的 2 -22范围内的完全可表示的结果，然后添加slo*slo并获得完全可表示的结果 ---s*s-foo精确计算。

当你除以 时s，你会产生额外的半 ulp 误差，这里最多为 2 -48，因为我们的误差已经很小了。

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

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

于 2013-07-11T19:51:23.427 回答