Я пытаюсь вычислить 32-битный корень с плавающей точкой IEEE-754 для различных входов, но для одного конкретного входа нижеприведенный алгоритм, основанный на методе Ньютона-Рафсона, выиграл 'не сходятся, мне интересно, что я могу сделать, чтобы решить эту проблему? Для платформы, которую я проектирую, у меня есть 32-битный сумматор / вычитатель с плавающей запятой, множитель и делитель.

Для ввода 0x7F7FFFFF (3.4028234663852886E38) алгоритм выигралt сходятся к правильному ответу 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)));