CORDIC help page gives inaccurate results

Question

Harry 2020-8-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/578949-cordic-help-page-gives-inaccurate-results

I am trying to use the code provided at this help page for computing square roots using CORDIC:

https://uk.mathworks.com/help/fixedpoint/examples/compute-square-root-using-cordic.html

I took the provided code, added the input initialization and output normalization (as described in the article) to get this function:

function s = cordic_sqrt(v, n)
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
% function s = cordic_sqrt(v, n)
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% Calculates the (approximate) square root of v over n CORDIC iterations.
%
% Taken from the Matlab help pages:
% mathworks.com/help/fixedpoint/examples/compute-square-root-using-cordic.html
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
% Set initial values
x = v + 0.25;
y = v - 0.25;
k = 4; % Used for the repeated (3*k + 1) iteration steps
for idx = 1:n
    xtmp = x * 2^(-idx);
    ytmp = y * 2^(-idx);
    if y < 0
        x = x + ytmp;
        y = y + xtmp;
    else
        x = x - ytmp;
        y = y - xtmp;
    end
    
    if idx==k
        xtmp = x * 2^(-idx);
        ytmp = y * 2^(-idx);
        if y < 0
            x = x + ytmp;
            y = y + xtmp;
        else
            x = x - ytmp;
            y = y - xtmp;
        end
        k = 3*k + 1;
    end
end
% Calculate normalization
An = prod(sqrt(1 - 2.^(-2*(1:n))));
% Normalize
s = x / An;

However, this gives very inaccurate results. For example:

cordic_sqrt(3/8, 30)
ans =
         0.611175220934138

(whereas MATLAB's sqrt function gives sqrt(3/8) = 0.612372435695794). I notice that if I comment out this part of the code:

%     if idx==k
%         xtmp = x * 2^(-idx);
%         ytmp = y * 2^(-idx);
%         if y < 0
%             x = x + ytmp;
%             y = y + xtmp;
%         else
%             x = x - ytmp;
%             y = y - xtmp;
%         end
%         k = 3*k + 1;
%     end

then I get accurate results (identical to MATLAB's sqrt function). So what is going on here?

I feel suspicious about the quality of the code provided in the article. I don't think the definition of the normalization term

is exactly correct (since it doesn't take into account these extra "if idx==k" iterations). However, this won't make any difference in my example because my number of iterations is so high (30).