CORDIC for square roots

Question

I have been looking at the CORDIC algorithm in hyperbolic rotation mode to find the square root of a variable. I am unsure what my initial variables should be (x0, y0, z0) to find the square root. I have read some papers, citing that to find the sqrt(a), initial values should be set to a+1,a-1,0 for x0,y0,and z0 respectively. Others says it should be a+0.25,a-0.25,0. I am very confused by this

Can anyone help?

double x = (64.0+1);
double y = (64.0-1);
double z = 0; 


double k = 3;
double n = 1;
while(n <= 20 ){

    double xn = pow(2.0,-1.0*n) * x;
    double yn = pow(2.0,-1.0*n) * y;

    if(y < 0){ 
        x = x + xn;
        y = y + yn;
        z = z - atanh(pow(2.0,-1.0*n));
    }
    else
    {
        x = x - xn;
        y = y - yn;
        z = z + atanh(pow(2.0,-1*n));

    }

    if(k > 0){
        k = k-1;
    }
    else{
        k = 3;
        if(y < 0){ 
            x = x + xn;
            y = y + yn;
            z = z - atanh(pow(2.0,-1.0*n));
        }
        else
        {
            x = x - xn;
            y = y - yn;
            z = z + atanh(pow(2.0,-1.0*n));

        }
    }
    n++;
    cout << "x: " << x << " y: " << y << " z: " << z << endl;
}

EDIT* Along with compensating for 3j+1 repeats, CORDIC requires to execute the loop twice in instances such as n = 4,13,40, ... I have updated my code to compensate for that, but it still does not work. I am using hyperbolic rotation in vectoring mode, which the variable d should be based on sign of y

EDIT* Turns out that CORDIC can fail when computing larger square root values, so you have to normalize the number you are trying to find the square root of to 0.5 to 2 range, then scale back answer back up.

Answer 1

You should normalize the number you are performing the square root on to [0.5,2) range, then scale back up accordingly.

Answer 2

initial values should be set to a+1,a-1,0 for x0,y0,and z0 respectively. Others says it should be a+0.25,a-0.25,0. I am very confused by this

The final result is sqrt((a+1)^2 - (a-1)^2) or sqrt((a+0.25)^2 - (a-0.25)^2). Either way the a^2 terms cancel and the constant terms cancel. The only difference is the first version returns sqrt(4a) or 2sqrt(a) and the second one returns sqrt(a) directly. I don't know the numerical reasons why one case or the other might be preferred.

Edit: Your bug is setting d based on y, it should be based on z.

Answer 3

There are several bugs in the code!

1. The variables are not exchanged.

     x = x + xn;
     y = y + yn;

should be:

     x = x + yn;
     y = y + xn;

2. The iteration 4 and 13 are not properly repeated. The following needs to be recomputed in those branches:

     xn = pow(2.0,-1.0*n) * y;
     yn = pow(2.0,-1.0*n) * x;

3. Scale correction should be done on the output (x/2*1.207497).

See here for ugly debugged, but working code: https://godbolt.org/z/ma4dCV

Note also, there is no need for the whole z channel (https://mathworks.com/help/fixedpoint/examples/compute-square-root-using-cordic.html).