Compute Square Root Using CORDIC

This example shows how to compute square root using a CORDIC kernel algorithm in MATLAB®. CORDIC-based algorithms are critical to many embedded applications, including motor controls, navigation, signal processing, and wireless communications.

CORDIC is an acronym for COordinate Rotation DIgital Computer. The Givens rotation-based CORDIC algorithm is one of the most hardware efficient algorithms because it only requires iterative shift-add operations [1,2]. The CORDIC algorithm eliminates the need for explicit multipliers and is suitable for calculating a variety of functions, such as sine, cosine, arcsine, arccosine, arctangent, vector magnitude, divide, square root, hyperbolic, and logarithmic functions.

The fixed-point CORDIC algorithm requires the following operations:

• 1 table lookup per iteration

• 2 shifts per iteration

• 3 additions per iteration

Note that for hyperbolic CORDIC-based algorithms, such as square root, certain iterations ($$i=4,13,40,121,\dots ,k,3k+1,\dots$$) are repeated to achieve result convergence [2].

CORDIC Kernel Algorithms Using Hyperbolic Computation Modes

You can use a CORDIC computing mode algorithm to calculate hyperbolic functions, such as hyperbolic trigonometric, square root, log, exp, etc.

CORDIC Equations in Hyperbolic Vectoring Mode

The hyperbolic vectoring mode is used for computing square root.

For the vectoring mode, the CORDIC equations are as follows:

$$\begin{array}{lcl}{x}_{i+1}& =& x_i+y_i\cdot d_i\cdot 2^{-i}\\ y_{i+1}& =& y_i+x_i\cdot d_i\cdot 2^{-i}\\ z_{i+1}& =& z_i-d_i\cdot \text{atanh}(2^{-i})\end{array}$$

where

$$i=\{1,2,3,4,4,5,\dots \}$$ is the sequence with repeated steps every $$3k+1$$iteration, and

$d_i=+1$ if $y_i<0$, and $-1$ otherwise.

This mode provides the following result as $N$ approaches $+\infty$:

$$\begin{array}{lcl}{x}_N& \approx & A_N\sqrt{x_0^2-y_0^2}\\ y_N& \approx & 0\\ z_N& \approx & z_0+\text{atanh}(y_0/x_0)\end{array}$$

where

$A_N={\prod }_{i=\{1,2,3,4,4,5,\dots \}}\sqrt{1-2^{-2i}}$.

Note that the sequence in the product is repeated on each of the $$3k+1$$ steps [2].

Typically $N$ is chosen to be a large-enough constant value. Thus, $A_N$ may be pre-computed.

Note also that for square root, only the $x_N$ result is used.

Implement a CORDIC Hyperbolic Vectoring Algorithm in MATLAB

A MATLAB code implementation example of the CORDIC hyperbolic vectoring kernel algorithm follows (for the case of scalar x, y, and z). This same code can be used for both fixed-point and floating-point data types.

CORDIC Hyperbolic Vectoring Kernel

k = 4; % Used for the repeated (3*k + 1) iteration steps
for idx = 1:N
    xtmp = bitsra(x, idx); % multiply by 2^(-idx)
    ytmp = bitsra(y, idx); % multiply by 2^(-idx)
    if y < 0
        x(:) = x + ytmp;
        y(:) = y + xtmp;
        z(:) = z - atanhLookupTable(idx);
    else
        x(:) = x - ytmp;
        y(:) = y - xtmp;
        z(:) = z + atanhLookupTable(idx);
    end
    if idx==k
        xtmp = bitsra(x, idx); % multiply by 2^(-idx)
        ytmp = bitsra(y, idx); % multiply by 2^(-idx)
        if y < 0
            x(:) = x + ytmp;
            y(:) = y + xtmp;
            z(:) = z - atanhLookupTable(idx);
        else
            x(:) = x - ytmp;
            y(:) = y - xtmp;
            z(:) = z + atanhLookupTable(idx);
        end
        k = 3*k + 1;
    end
end % idx loop

The hyperbolic inverse tangent lookup table atanhLookupTable is defined as follows, which in code generation is pre-computed and will be a constant.

atanhLookupTable = cast(atanh(2.^-(1:N)),'like',z);

Compute Square Root Using the CORDIC Hyperbolic Vectoring Kernel

The judicious choice of initial values allows the CORDIC kernel hyperbolic vectoring mode algorithm to compute square root.

First, the following initialization steps are performed:

• $x_0$ is set to $v+0.25$.

• $y_0$ is set to $v-0.25$.

After $N$ iterations, these initial values lead to the following output as $N$ approaches $+\infty$:

$$x_N\approx A_N\sqrt{(v+0.25{)}^2-(v-0.25{)}^2}$$

This may be further simplified as follows:

$$x_N\approx A_N\sqrt{v}$$

where $A_N$ is the CORDIC gain as defined above.

Note that for square root, $z$ and atanhLookupTable have no impact on the result. Hence, $z$ and atanhLookupTable are not used.

MATLAB Implementation of a CORDIC Square Root Kernel

A MATLAB code implementation example of the CORDIC square root kernel algorithm follows (for the case of scalar x and y). This same code can be used for both fixed-point and floating-point data types.

CORDIC Square Root Kernel

k = 4; % Used for the repeated (3*k + 1) iteration steps
for idx = 1:N
    xtmp = bitsra(x, idx); % multiply by 2^(-idx)
    ytmp = bitsra(y, idx); % multiply by 2^(-idx)
    if y < 0
        x(:) = x + ytmp;
        y(:) = y + xtmp;
    else
        x(:) = x - ytmp;
        y(:) = y - xtmp;
    end
    if idx==k
         xtmp = bitsra(x, idx); % multiply by 2^(-idx)
         ytmp = bitsra(y, idx); % multiply by 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 % idx loop

This code is identical to the CORDIC hyperbolic vectoring kernel implementation above, except that z and atanhLookupTable are not used. This is a cost savings of 1 table lookup and 1 addition per iteration.

Example

Use the cordicsqrt function to compute the approximate square root of v_fix using ten CORDIC kernel iterations.

step  = 2^-7;
v_fix = fi(0.5:step:(2-step), 1, 20); % fixed-point inputs in range [.5, 2)
niter = 10; % number of CORDIC iterations
x_sqr = cordicsqrt(v_fix, niter);

Get the real-world value (RWV) of the CORDIC outputs for comparison and plot the error between the MATLAB reference and CORDIC square root values.

x_cdc = double(x_sqr); % CORDIC results (scaled by An_hp)
v_ref = double(v_fix); % Reference floating-point input values
x_ref = sqrt(v_ref);   % MATLAB reference floating-point results
figure;
subplot(211);
plot(v_ref, x_cdc, 'r.', v_ref, x_ref, 'b-');
legend('CORDIC', 'Reference', 'Location', 'SouthEast');
title('CORDIC Square Root (In-Range) and MATLAB Reference Results');
subplot(212);
absErr = abs(x_ref - x_cdc);
plot(v_ref, absErr);
title('Absolute Error (vs. MATLAB SQRT Reference Results)');

Overcoming Algorithm Input Range Limitations

Many square root algorithms normalize the input value, $v$, to within the range of [0.5, 2). This pre-processing is typically done using a fixed word length normalization, and can be used to support small as well as large input value ranges.

The CORDIC-based square root algorithm implementation is particularly sensitive to inputs outside of this range. The cordicsqrt function overcomes this algorithm range limitation through a normalization approach based on the following mathematical relationships:

$v=u\cdot 2^n$, for some $0.5<=u<2$ and some even integer $n$.

Thus:

$$\sqrt{v}=\sqrt{u}\cdot 2^{n/2}$$

In the cordicsqrt function, the values for $u$ and $n$, described above, are found during normalization of the input $v$. $n$ is the number of leading zero most-significant bits (MSBs) in the binary representation of the input $v$. These values are found through a series of bit-wise logic and shifts. Note that because $n$ must be even, if the number of leading zero MSBs is odd, one additional bit shift is made to make $n$ even. The resulting value after these shifts is the value $0.5<=u<2$.

$u$ becomes the input to the CORDIC-based square root kernel, where an approximation to $\sqrt{u}$ is calculated. The result is then scaled by $2^{n/2}$ so that it is back in the correct output range. This is achieved through a simple bit shift by $n/2$ bits. The (left or right) shift direction depends on the sign of $n$.

Example

Compute the square root of 10-bit fixed-point input data with a small non-negative range using CORDIC. Compare the CORDIC-based algorithm results to the floating-point MATLAB reference results over the same input range.

step     = 2^-8;
u_ref    = 0:step:(0.5-step); % Input array (small range of values)
u_in_arb = fi(u_ref,0,10);    % 10-bit unsigned fixed-point input data values
u_len    = numel(u_ref);
sqrt_ref = sqrt(double(u_in_arb)); % MATLAB sqrt reference results
niter    = 10;
results  = zeros(u_len, 2);
results(:,2) = sqrt_ref(:);

Compute the equivalent real-world value (RWV) result. Plot the RWV of CORDIC and MATLAB reference results.

x_out = cordicsqrt(u_in_arb, niter);
results(:,1) = double(x_out);
figure;
subplot(211);
plot(u_ref, results(:,1), 'r.', u_ref, results(:,2), 'b-');
legend('CORDIC', 'Reference', 'Location', 'SouthEast');
title('CORDIC Square Root (Small Input Range) and MATLAB Reference Results');
axis([0 0.5 0 0.75]);
subplot(212);
absErr = abs(results(:,2) - results(:,1));
plot(u_ref, absErr);
title('Absolute Error (vs. MATLAB SQRT Reference Results)');

Example

Compute the square root of 16-bit fixed-point input data with a large positive range using CORDIC. Compare the CORDIC-based algorithm results to the floating-point MATLAB reference results over the same input range.

u_ref    = 0:5:2500;       % Input array (larger range of values)
u_in_arb = fi(u_ref,0,16); % 16-bit unsigned fixed-point input data values
u_len    = numel(u_ref);
sqrt_ref = sqrt(double(u_in_arb)); % MATLAB sqrt reference results
niter    = 16;
results  = zeros(u_len, 2);
results(:,2) = sqrt_ref(:);

Compute the equivalent real-world value (RWV) result. Plot the RWV of CORDIC and MATLAB reference results.

x_out = cordicsqrt(u_in_arb, niter);
results(:,1) = double(x_out);
figure;
subplot(211);
plot(u_ref, results(:,1), 'r.', u_ref, results(:,2), 'b-');
legend('CORDIC', 'Reference', 'Location', 'SouthEast');
title('CORDIC Square Root (Large Input Range) and MATLAB Reference Results');
axis([0 2500 0 55]);
subplot(212);
absErr = abs(results(:,2) - results(:,1));
plot(u_ref, absErr);
title('Absolute Error (vs. MATLAB SQRT Reference Results)');

References

1. Jack E. Volder, "The CORDIC Trigonometric Computing Technique," IRE Transactions on Electronic Computers, Volume EC-8, September 1959, pp. 330-334.

2. J.S. Walther, "A Unified Algorithm for Elementary Functions," Conference Proceedings, Spring Joint Computer Conference, May 1971, pp. 379-385.