CORDIC

CORDIC is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Introduction

CORDIC is the name of an algorithm for calculating trigonometric, logarithmic and hyperbolic functions, named after its first application on an airborne computer (COordinate Rotation DIgital Computer) in 1959. Unlike a Taylor expansion or polynomial approximation, it converges rapidly on machines with low computing and memory capacities: to calculate a tangent with 10 significant digits, it requires only 6 floating-point constants, and only additions, subtractions and digit shifts in its iterative part.

It is valid for angle values between 0 and π/2 only, but whatever the value of an angle, the calculation of its tangent can always be reduced to that of an angle between 0 and π/2, using trigonometric identities. Similarly, once you know the tangent, you can easily calculate the sine or cosine.

Pseudo code

constant θ[n] = arctan 10^(-n) // or simply 10^(-n) depending on floating point precision 
constant epsilon = 10^-12

function tan(alpha)            // 0 < alpha <= π/2 
  x = 1 ; y = 0 ; k = 0
  while epsilon < alpha
    while alpha < θ[k] 
       k++
    end loop
    alpha -= θ[k]
    x2 = x - 10^(-k)*y
    y2 = y + 10^(-k)*x
    x = x2 ; y = y2
  end loop
  return (y/x)
end function

Task

• Implement the CORDIC algorithm, using only the 4 arithmetic operations and right shifts in the main loop if possible.
• Use your implementation to calculate the cosine of the following angles, expressed in radians: -9, 0, 1.5 and 6

See also

CORDIC on Wikipedia

ALGOL 68

Based on the pseudo code, method to calculate cos (and sin) from the XPL0 sample.

BEGIN # implement sin, cos and tan using the CORDIC algorithm                #

    REAL pi by 2 = pi / 2;

    # pre-computed table of arctan(10^-n) values:                            #
    []REAL theta = ( 7.85398163397448e-01, 9.96686524911620e-02
                   , 9.99966668666524e-03, 9.99999666666867e-04
                   , 9.99999996666667e-05, 9.99999999966667e-06
                   , 9.99999999999667e-07, 9.99999999999997e-08
                   , 1.00000000000000e-08, 1.00000000000000e-09
                   , 1.00000000000000e-10, 1.00000000000000e-11
                   , 1.00000000000000e-12, 1.00000000000000e-13
                   , 1.00000000000000e-14, 1.00000000000000e-15
                   , 1.00000000000000e-16
                   );
    REAL epsilon = 1e-16;

    # mode to hold the values returned from the CORDIC procedure             #
    MODE CORDICV = STRUCT( REAL y, x );

    # CORDIC algorithm, finds "y" and "x" for alpha radians                  #
    # signs indicates the sign of the results:                               #
    #     signs[ 1 ] = sign for -π : -π/2, signs[ 1 ] = sign for -π/2 : 0    #
    #     signs[ 3 ] = sign for  0 :  π/2, signs[ 4 ] = sign for  π/2 : π    #
    # sign both = TRUE => sign applied to both y and x, FALSE => sign y only #
    PROC cordic = ( REAL alpha in, []INT signs, BOOL sign both )CORDICV:
         BEGIN
            REAL alpha     := alpha in;            # ensure -π <= alpha <= π #
            BOOL flip sign := FALSE;
            WHILE alpha < - pi DO alpha +:= pi; flip sign := NOT flip sign OD;
            WHILE alpha >   pi DO alpha -:= pi; flip sign := NOT flip sign OD;
            INT    sign;
            IF     alpha < - pi by 2 THEN
                alpha +:= pi;
                sign := signs[ 1 ]
            ELIF   alpha < 0         THEN
                alpha := - alpha;
                sign := signs[ 2 ]
            ELIF   alpha < pi by 2   THEN
                sign := signs[ 3 ]
            ELSE # alpha < pi #
                alpha := pi - alpha;
                sign := signs[ 4 ]
            FI;
            IF flip sign AND sign both THEN sign := -sign FI;
            REAL x := 1, y := 0;
            INT  k := 1;
            REAL ten to minus k := 1;     # NB: ten to minus k = 10.0^-(k-1) #
            WHILE epsilon < alpha DO
                WHILE alpha < theta[ k ] DO
                    k              +:=  1;
                    ten to minus k /:= 10
                OD;
                alpha -:= theta[ k ];
                REAL x2 = x - ten to minus k * y;
                REAL y2 = y + ten to minus k * x;
                x      := x2;
                y      := y2
            OD;
            CORDICV( sign * y, IF sign both THEN sign * x ELSE x FI )
         END # cordic # ;

    # find sin(alpha) using the CORDIC algorithm. alpha in radians           #
    PROC c cos = ( REAL alpha )REAL:
         BEGIN
            CORDICV c = cordic( alpha, ( -1, 1, 1, -1 ), TRUE );
            x OF c / sqrt( ( x OF c * x OF c ) + ( y OF c * y OF c ) )
         END # c cos # ;

    # find sin(alpha) using the CORDIC algorithm. alpha in radians           #
    PROC c sin = ( REAL alpha )REAL:
         BEGIN
            CORDICV c = cordic( alpha, ( -1,-1, 1, 1 ), TRUE );
            y OF c / sqrt( ( x OF c * x OF c ) + ( y OF c * y OF c ) )
         END # c cos # ;

    # find tan(alpha) using the CORDIC algorithm. alpha in radians           #
    PROC c tan = ( REAL alpha )REAL:
         BEGIN
             CORDICV c = cordic( alpha, (  1, -1, 1,-1 ), FALSE );
             IF x OF c = 0 THEN max real ELSE y OF c / x OF c FI
         END # c tan # ;

    PROC show cordic = ( REAL angle )VOID:
         BEGIN
            REAL cosine  = cos( angle ), cordic cosine = c cos( angle );
            REAL sine    = sin( angle ), cordic sine   = c sin( angle );
            REAL tangent = tan( angle ), cordic tan    = c tan( angle );
            REAL c diff  = ABS ( cordic cosine - cosine  );
            REAL s diff  = ABS ( cordic sine   - sine    );
            REAL t diff  = ABS ( cordic tan    - tangent );
            print( ( fixed(         angle,  -8, 4    ), ":  "
                   , fixed( cordic cosine,  -9, 6    ), " "
                   , fixed(        cosine,  -9, 6    ), " "
                   , float(        c diff, -14, 8, 2 ), " | "
                   , fixed(   cordic sine,  -9, 6    ), " "
                   , fixed(          sine,  -9, 6    ), " "
                   , float(        s diff, -14, 8, 2 ), " | "
                   , fixed(    cordic tan,  -9, 6    ), " "
                   , fixed(       tangent,  -9, 6    ), " "
                   , float(        t diff, -14, 8, 2 ), newline
                   )
                 )
         END # show cordic # ;

    []REAL tests = ( -9, 0, 1.5, 6 );
    print( ( "    angle cordic cos       cos     difference" ) );
    print( (         "  cordic sin       sin     difference" ) );
    print( (         "  cordic tan       tan     difference" ) );
    print( ( newline ) );
    FOR i FROM LWB tests TO UPB tests DO
        show cordic( tests[ i ] )
    OD;
    print( ( "    angle cordic cos       cos     difference" ) );
    print( (         "  cordic sin       sin     difference" ) );
    print( (         "  cordic tan       tan     difference" ) );
    print( ( newline ) );
    FOR i FROM 7810 BY 9 TO 7898 DO
        show cordic( i / 10 000 )
    OD

END

Output:

    angle cordic cos       cos     difference  cordic sin       sin     difference  cordic tan       tan     difference
 -9.0000:  -0.911130 -0.911130  1.1102230e-16 | -0.412118 -0.412118  0.00000000e+0 |  0.452316  0.452316  0.00000000e+0
  0.0000:   1.000000  1.000000  0.00000000e+0 |  0.000000  0.000000  0.00000000e+0 |  0.000000  0.000000  0.00000000e+0
  1.5000:   0.070737  0.070737  7.4940054e-16 |  0.997495  0.997495  0.00000000e+0 | 14.101420 14.101420  1.5099033e-13
  6.0000:   0.960170  0.960170  1.1102230e-16 | -0.279415 -0.279415  5.5511151e-16 | -0.291006 -0.291006  6.1062266e-16
    angle cordic cos       cos     difference  cordic sin       sin     difference  cordic tan       tan     difference
  0.7810:   0.710210  0.710210  3.3306691e-16 |  0.703990  0.703990  2.2204460e-16 |  0.991242  0.991242  6.6613381e-16
  0.7819:   0.709576  0.709576  2.2204460e-16 |  0.704629  0.704629  1.1102230e-16 |  0.993028  0.993028  4.4408921e-16
  0.7828:   0.708942  0.708942  3.3306691e-16 |  0.705267  0.705267  4.4408921e-16 |  0.994817  0.994817  9.9920072e-16
  0.7837:   0.708307  0.708307  6.6613381e-16 |  0.705905  0.705905  5.5511151e-16 |  0.996609  0.996609  1.6653345e-15
  0.7846:   0.707671  0.707671  5.5511151e-16 |  0.706542  0.706542  3.3306691e-16 |  0.998405  0.998405  1.3322676e-15
  0.7855:   0.707035  0.707035  1.1102230e-16 |  0.707179  0.707179  1.1102230e-16 |  1.000204  1.000204  2.2204460e-16
  0.7864:   0.706398  0.706398  1.1102230e-16 |  0.707815  0.707815  2.2204460e-16 |  1.002006  1.002006  4.4408921e-16
  0.7873:   0.705761  0.705761  1.1102230e-16 |  0.708450  0.708450  1.1102230e-16 |  1.003811  1.003811  2.2204460e-16
  0.7882:   0.705123  0.705123  3.3306691e-16 |  0.709085  0.709085  5.5511151e-16 |  1.005619  1.005619  1.3322676e-15
  0.7891:   0.704484  0.704484  5.5511151e-16 |  0.709720  0.709720  4.4408921e-16 |  1.007431  1.007431  1.3322676e-15

C

Translation of: Wren

If the M_PI constant is not defined in your math.h, then you can use 4 * atan(1) instead.

#include <stdio.h>
#include <math.h>

/* The following are pre-computed to avoid using atan and sqrt functions. */
const double angles[] = {
    0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
    0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
    0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119,
    0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812,
    0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863,
    0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929,
    0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058
};

const double kvalues[] = {
    0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775,
    0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889,
    0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894,
    0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314,
    0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925,
    0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888
};

double radians(double degrees) { return degrees * M_PI / 180.0; }

void Cordic(double alpha, int n, double *c_cos, double *c_sin) {
    int i, ix, sigma;
    double kn, x, y, atn, t, theta = 0.0, pow2 = 1.0;
    int newsgn = (int)floor(alpha / (2.0 * M_PI)) % 2 == 1 ? 1 : -1;
    if (alpha < -M_PI/2.0 || alpha > M_PI/2.0) {
        if (alpha < 0) {
            Cordic(alpha + M_PI, n, &x, &y);
        } else {
            Cordic(alpha - M_PI, n, &x, &y);
        }
        *c_cos = x * newsgn;
        *c_sin = y * newsgn;
        return;
    }
    ix = n - 1;
    if (ix > 23) ix = 23;
    kn = kvalues[ix];
    x = 1;
    y = 0;
    for (i = 0; i < n; ++i) {
        atn = angles[i];
        sigma = (theta < alpha) ? 1 : -1;
        theta += sigma * atn;
        t = x;
        x -= sigma * y * pow2;
        y += sigma * t * pow2;
        pow2 /= 2.0;
    }
    *c_cos = x * kn;
    *c_sin = y * kn;
}

int main() {
    int i, th;
    double thr, c_cos, c_sin;
    double angles[] = {-9.0, 0.0, 1.5, 6.0 };
    char *f;
    printf("  x       sin(x)     diff. sine     cos(x)    diff. cosine\n");
    f =  "%+03d.0°  %+.8f (%+.8f) %+.8f (%+.8f)\n";
    for (th = -90; th <= +90; th += 15) {
        thr = radians(th);
        Cordic(thr, 24, &c_cos, &c_sin);
        printf(f, th, c_sin, c_sin - sin(thr), c_cos, c_cos - cos(thr));
    }
    printf("\nx(rads)   sin(x)     diff. sine     cos(x)    diff. cosine\n");
    f = "%+4.1f    %+.8f (%+.8f) %+.8f (%+.8f)\n";
    for (i = 0; i < 4; ++i) {
        thr = angles[i];
        Cordic(thr, 24, &c_cos, &c_sin);
        printf(f, thr, c_sin, c_sin - sin(thr), c_cos, c_cos - cos(thr));
    }
    return 0;
}

Output:

  x       sin(x)     diff. sine     cos(x)    diff. cosine
-90.0°  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
-75.0°  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
-60.0°  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
-45.0°  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
-30.0°  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
-15.0°  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
+00.0°  +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+15.0°  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
+30.0°  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
+45.0°  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
+60.0°  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
+75.0°  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
+90.0°  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

x(rads)   sin(x)     diff. sine     cos(x)    diff. cosine
-9.0    -0.41211842 (+0.00000006) -0.91113029 (-0.00000003)
+0.0    +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+1.5    +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
+6.0    -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)

FreeBASIC

Translation of: Julia

#define min(a, b)  iif((a) < (b), (a), (b))
#define floor(x)  ((x*2.0-0.5) Shr 1)
#define pi  4 * Atn(1)

Sub CORDIC(alfa As Double, iteration As Integer, Byref cosAlfa As Double, Byref sinAlfa As Double)
    ' Fix for the Wikipedia's MATLAB code bug in sin (sometimes flips sign) when |?| > 2p
    Dim As Integer newsgn = Iif(Int(floor(alfa / (2 * pi))) Mod 2 = 1, 1, -1)
    If alfa < -Atn(1) * 2 Then
        CORDIC(alfa + pi, iteration, cosAlfa, sinAlfa)
        cosAlfa *= newsgn
        sinAlfa *= newsgn
        Exit Sub
    End If
    If alfa > Atn(1) * 2 Then
        CORDIC(alfa - pi, iteration, cosAlfa, sinAlfa)
        cosAlfa *= newsgn
        sinAlfa *= newsgn
        Exit Sub
    End If
    
    ' Initialization of tables of constants used by CORDIC
    ' need a table of arctangents of negative powers of two, in radians:
    ' angles = atan(2.^-(0:27));
    Dim As Double angles(27) = { _
    0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676, _
    0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010, _
    0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119, _
    0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812, _
    0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863, _
    0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929, _
    0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058 }
    
    ' and a table of products of reciprocal lengths of vectors (1, 2^-2j):
    ' Kvalores = cumprod(1./sqrt(1 + 1j*2.^(-(0:23))))
    Dim As Double Kvalues(23) = { _
    0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775, _
    0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889, _
    0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894, _
    0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314, _
    0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925, _
    0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888 }
    
    Dim As Double Kn = Kvalues(min(iteration, Ubound(Kvalues)))
    Dim As Double v(1) = {1, 0} ' start with 2-vector cosine and sine of zero
    Dim As Double poweroftwo = 1
    Dim As Double angle = angles(0)
    
    ' Iterations
    For j As Integer = 0 To iteration - 1
        Dim As Integer sigma = Iif(alfa < 0, -1, 1)
        Dim As Double factor = sigma * poweroftwo
        ' Note the matrix multiplication can be done using scaling by powers of two and addition subtraction
        Dim As Double R(1 To 2, 1 To 2)
        R(1, 1) = 1 :        R(1, 2) = -factor
        R(2, 1) = factor :   R(2, 2) = 1
        Dim As Double temp(1)
        temp(0) = R(1, 1) * v(0) + R(1, 2) * v(1)
        temp(1) = R(2, 1) * v(0) + R(2, 2) * v(1)
        v(0) = temp(0)
        v(1) = temp(1)
        alfa -= sigma * angle ' update the remaining angle
        poweroftwo /= 2
        ' update the angle from table, or eventually by just dividing by two
        If j + 1 >= Ubound(angles) Then
            angle /= 2
        Else
            angle = angles(j + 1)
        End If
    Next
    
    ' Adjust length of output vector to be [cos(alfa), sin(alfa)]
    cosAlfa = v(0) * Kn
    sinAlfa = v(1) * Kn
End Sub

Print "x(deg)    sin(x)       diff. sine   cos(x)       diff. cosine "
Dim As Double cosTheta, sinTheta, theta ', thetaR
For theta = -90 To 90 Step 15
    CORDIC(theta * Atn(1) / 45, 24, cosTheta, sinTheta)
    Print Using "+##.#" & Chr(248) & "  +#.######## (+#.########) +#.######## (+#.########)"; theta; sinTheta; sinTheta - Sin(theta * Atn(1) / 45); cosTheta; cosTheta - Cos(theta * Atn(1) / 45)
Next

Print !"\nx(rad)    sin(x)       diff. sine   cos(x)       diff. cosine "
Dim As Double thetaR(1 To 4) = {-9, 0, 1.5, 6}
For r As Integer = 1 To 4 '-9 To 6 Step 1.5
    CORDIC(thetaR(r), 24, cosTheta, sinTheta)
    Print Using " +#.#   +#.######## (+#.########) +#.######## (+#.########)"; thetaR(r); sinTheta; sinTheta - Sin(thetaR(r)); cosTheta; cosTheta - Cos(thetaR(r))
Next

Sleep

Output:

x(deg)    sin(x)       diff. sine   cos(x)       diff. cosine
-90.0°  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
-75.0°  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
-60.0°  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
-45.0°  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
-30.0°  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
-15.0°  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
 +0.0°  -0.00000007 (-0.00000007) +1.00000000 (-0.00000000)
+15.0°  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
+30.0°  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
+45.0°  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
+60.0°  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
+75.0°  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
+90.0°  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

x(rad)    sin(x)       diff. sine   cos(x)       diff. cosine
 -9.0   -0.41211842 (+0.00000006) -0.91113029 (-0.00000003)
 +0.0   -0.00000007 (-0.00000007) +1.00000000 (-0.00000000)
 +1.5   +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
 +6.0   -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)

Imp77

! This is effectively a fixed-point implementation of CORDIC sin and cos which
! checks results against the floating-point library used by C. There are a few
! things we could do to improve accuracy - but this is good enough for a quick
! hack that does not use any external math library code, other than the afore-
! mentioned consistency checks.

! The program is based on Michael Bertrand's CORDIC article which has appeared
! in multiple places including Dr Dobbs magazine, a copy of which can be found
! at http://www.pennelynn.com/Documents/CUJ/HTML/92HTML/1992024C.HTM

! Note that in creating this implementation, we hit a bug that manifested at
! the the extremes of the sin and cos functions (where the values are close
! to 0 or 1) and which I have patched with an approximation since the functions
! are close to linear at those points.  However the hack to correct that error
! does make this code a little larger than the C version from which it was
! derived.

! I have thrown in a local implementation of atan and sqrt for good measure :-)

! Imp77 does not have unsigned integer data types; however since the CORDIC code
! effectively uses a 1.14 fixed point implementation to implement angles (with
! the original C sources being compatible with a 16-bit int environment), and
! the current Imp77 compiler being a 32 bit implementation, unsigned integers
! have not required.

%begin

%const %integer debug = 0 { 0 : 4 }

!%external %real %fn %spec atan %alias "atanf" (%real angle)  ;! borrow C's math library for use in demo.
%real %fn atan(%longreal x1)
  %const %long %real r1= 1.2595802263029547 @  1
  %const %long %real r2=-8.6186317517509520 @  1
  %const %long %real r3=-1.2766919133361079 @  0
  %const %long %real r4=-8.3921038065840512 @ -2
  %const %long %real s1= 2.7096164294378656 @  1
  %const %long %real s2= 6.5581320451487386 @  0
  %const %long %real s3= 2.1441643116703661 @  0
  %const %long %real s4= 1.2676256708212610 @  0
  %const %long %real rt3= 1.7320508075688772 @  0
  %const %long %real piby6= 5.2359877559829887 @ -1
  %const %long %real piby2m1= 5.7079632679489661 @ -1
  %const %long %real rt3m1=7.3205080756887728 @ -1
  %const %long %real tanpiby12= 2.6794919243112271 @ -1
  %long %real xx1,xsq,constant=0
  %integer sign=0,inv=0
  
  %if x1<0 %then sign=1 %and xx1=-x1 %else xx1=x1
  %if xx1>1.0 %then xx1=1.0/xx1 %and inv=1
  %if xx1>tanpiby12 %then xx1=(rt3m1*xx1-1.0+xx1)/(xx1+rt3) %and constant=piby6
  xsq=xx1*xx1
  xx1 = xx1*(r1/(xsq+s1+r2/(xsq+s2+r3/(xsq+s3+r4/(xsq+s4)))))
  xx1 = xx1 + constant
  %if inv=1 %then xx1=1.0-xx1+piby2m1
  %if sign=1 %then xx1=-xx1
  %result=xx1 { shorten from %long %real down to %real when returning result }
%end

!%external %real %fn %spec sqrt %alias "sqrtf" (%real angle)
%real %fn sqrt(%long %real X)
  %integer %fn approx integer sqrt(%integer X)
    %integer k = 1, result = 1
    k = k+1 %and result = k*k %while result <= X
    %result = k-1 {note 'off by 1' bug in original C version due to returning k--}
  %end
  %integer I
  %long %real prev A, A, epsilon = X / 100  {accurate enough due to .14 fixed point}
  %if X = 2.0 %then %result = 1.41421356237 %else I = approx integer sqrt(INT PT(X))
  %result = I %if X = I*I { exact root }
  ! Newton's method: A{n} = (A{n-1} + X/A{n-1})/2
  A = I; prev A = A %and A = (A + X/A)/2.0 %until prev A - epsilon <= A %and A <= prev A + epsilon
  %result = A
%end

%real PI = atan(1) * 4

%const %integer K64 = 1<<16, K32 = K64>>1, K16 = K32>>1, K8  = K16>>1

%routine phex(%integer x)
  %const %string (16) hex = "0123456789ABCDEF"
  %integer i
  print symbol(charno(hex, (x >> (i<<2))&15+1)) %for i = 7, -1, 0
%end

%routine SinCos({unsigned} %integer theta, %integer %name sin, %integer %name cos)
  %const %integer QUAD1 = 3, QUAD2 = 2, QUAD3 = 0, QUAD4 = 1   { bit 1 for COS x vs -x,  bit 2 for SIN x vs -x }
  %const %integer NBITS = 14                                   { NBITS is number of bits for CORDIC units. }

  {unsigned} %const %integer CordicBase = 1 << NBITS                      { base for CORDIC units }
  {unsigned} %const %integer HalfBase = CordicBase >> 1                   { CordicBase / 2 }
  {unsigned} %const %integer Quad2Boundary = CordicBase << 1              { 2 * CordicBase }
  {unsigned} %const %integer Quad3Boundary = CordicBase + Quad2Boundary   { 3 * CordicBase }

  ! these owns are initialised on the first call to SinCos:
  %own %integer %array ArcTan(0:NBITS-1)                       { angles for rotation }
  %own %integer xInit                                          { initial x projection }
  
  %integer xcos,xsin  { used to patch a problem that happens near the zeroes and extremes of the results }
  
  !  USE : Calc sin and cos with integer CORDIC routine.
  !  IN  : theta = incoming angle (in CORDIC angle units)
  !  OUT : sin (in CORDIC fixed point units)
  !        cos (in CORDIC fixed point units)
  !  NOTE: The incoming angle theta is in CORDIC angle
  !        units, which subdivide the circle into 64K parts,
  !        with 0 deg = 0, 90 deg = 16384 (CordicBase), 180 deg
  !        = 32768, 270 deg = 49152, etc. The calculated sine
  !        and cosine are in CORDIC fixed point units - an int
  !        considered as a fraction of 16384 (CordicBase).

  %integer quadrant  { quadrant of incoming angle }
  %integer z         { incoming angle moved to 1st quad }
  %integer i         { to index rotations - one per bit }
  %integer x, y      { projections onto axes }
  %integer xl, yl    { projections of rotated vector }

  %own %integer initialised = 0
  %if initialised = 0 %start

    !   USE : Load globals used by SinCos().
    !   OUT : Loads globals used in SinCos() :
    !    CordicBase    = base for CORDIC units
    !    HalfBase      = Cordicbase / 2
    !    Quad2Boundary = 2 * CordicBase
    !    Quad3Boundary = 3 * CordicBase
    !    ArcTan()      = the arctans of 1/(2^i)
    !    xInit         = initial value for x projection
    !   NOTE: Must be called once only to initialize before
    !    calling SinCos(). xInit is sufficiently less than
    !    CordicBase to exactly compensate for the expansion
    !    in each rotation.
  
    {%integer i}      { to index ArcTan() }
    %real f           { to calc initial x projection }
    %integer powr     { powers of 2 up to 2^(2*(NBITS-1)) }
  
    { ArcTans are diminishingly small angles. }
    powr = 1
    %for i = 0, 1, NBITS-1 %cycle
    
      ! We could alternatively have created ArcTan as an array of constants at compile time.
      ArcTan(i) = INT(atan(1.0/powr)/(PI/2.0)*CordicBase)
      %if debug > 1 %start
        print string("Setting ArcTan("); write(i,0); print string(") to "); phex(ArcTan(i)); newline
      %finish
      powr = powr << 1
      ! xInit is initial value of x projection to compensate for expansions. f = 1/sqrt(2/1 * 5/4 * ...
      ! Normalize as an NBITS binary fraction (multiply by CordicBase) and store in xInit.
      ! Get f = 0.607253 and xInit = 9949 = 0x26DD for NBITS = 14.
    %repeat
    f = 1.0

    powr = 1
    %for i = 0, 1, NBITS-1 %cycle
      f = (f * (powr + 1)) / powr
      powr = powr << 2
    %repeat

    f = 1.0/sqrt(f);              { sqrt is expensive but only done once during initialisation }
                                  { Note could have implemented f = rsqrt(f) instead and been faster }
                                  
    xInit = INT(CordicBase * f)   { rounding type conversion to fixed point }
    %if debug # 0 %start
      print string("xInit = "); print(CordicBase * f, 12)
      print string(" = "); write(xInit, 0); newlines(2)
    %finish
    
    initialised = 1
  %finish  { initialisation }

  ! Determine quadrant of incoming angle, translate to 1st quadrant.
  ! Note use of previously calculated values CordicBase, etc. for speed.

  theta = theta & (K64-1)  { reduce to a single circle }
  %if       theta < CordicBase    %then quadrant = QUAD1 %and z = theta                 %c
  %else %if theta < Quad2Boundary %then quadrant = QUAD2 %and z = Quad2Boundary - theta %c
  %else %if theta < Quad3Boundary %then quadrant = QUAD3 %and z = theta - Quad2Boundary %c
  %else                                 quadrant = QUAD4 %and z = -theta

  ! Initialize projections; Negate z, so same rotations taking angle z to 0 will take (x, y) = (xInit, 0) to (cos, sin).
  z = z & (K32-1) { this is a bit dodgy on my part as there *are* cases where z is negative above }
  x = xInit; y = 0; z = -z

  %for i = 0, 1, NBITS-1 %cycle                                                    { Rotate NBITS times. }
    %if z < 0 %then z = z+ArcTan(i) %and yl = y + (x >> i) %and xl = x - (y >> i)  { Counter-clockwise rotation. } %c
    %else z = z-ArcTan(i) %and yl = y - (x >> i) %and xl = x + (y >> i)            { Clockwise rotation. }
    x = xl; y = yl                                                                 { Put new projections into (x,y) for next go. }
  %repeat

  %if quadrant&1 # 0 %then cos = x %else cos = -x  { Attach signs depending on quadrant. }
  %if quadrant&4 # 0 %then sin = y %else sin = -y

  xsin=sin {debug}
  sin = -sin %if Theta < K32
  sin = sin & (K64-1); sin = sin!(\(K16-1)) %if sin > K16

  xcos=cos {debug}
  cos = cos & (K64-1); cos = cos!(\(K16-1)) %if cos > K16

  %if debug > 1 %start
    print string("                  Theta ="); write(theta*1608//1024,8)
    print string("; SIN(Theta) = "); phex(xsin);
    print string("; COS(Theta) = "); phex(xcos);
    newline
  %finish

  ! There is a bug in this implementation that fails near multiples of 90 degrees -
  ! this test catches the broken cases and uses an approximation to repair the results :-(
  ! I will remove these hacks later, if I ever find and fix the bug properly!
  ! My suspicion is of overflow of a partial 16x16 multiplication somewhere, or maybe a '<<'.
  %if xsin&16_FFFFC000 # 16_FFFFC000 %and xsin&16_FFFFC000 # 0 %then %start
    %if                  theta < 16_0400 %start
      sin = theta*1608//1024
    %else %if 16_3C00 <= theta < 16_4000
      sin = K16 - (16_4000-theta)//25
    %else %if 16_4000 <= theta < 16_4400
      sin = K16 - (theta-16_4000)//25
    %else %if 16_7C00 <= theta < 16_8400
      sin = (16_8000-theta)*1608//1024
    %else %if 16_BC00 <= theta < 16_C000
      sin = -K16 + (16_C000-theta)//25
    %else %if 16_C000 <= theta < 16_C400
      sin = -K16 + (theta-16_C000)//25
    %else %if 16_FC00 <= theta
      sin = -(16_10000-theta)*1608//1024
    %else
      print string("Oops, missed a case: sin "); phex(theta); newline; %stop
    %finish
  %finish
  
  %if xcos&16_FFFFC000 # 16_FFFFC000 %and xcos&16_FFFFC000 # 0 %then %start
    %if                  theta < 16_0400 %start
      cos = K16 - (theta)//25
    %else %if 16_3C00 <= theta < 16_4400
      cos = (16_4000-theta)*1608//1024
    %else %if 16_7C00 <= theta < 16_8000
      cos = -K16 + (16_8000-theta)//25
    %else %if 16_8000 <= theta < 16_8400
      cos = -K16 + (theta-16_8000)//25
    %else %if 16_BC00 <= theta < 16_C400
      cos = (theta-16_C000)*1608//1024
    %else %if 16_FC00 <= theta
      cos = K16 - (16_10000-theta)//25
    %else
      print string("Oops, missed a case: cos "); phex(theta); newline; %stop
    %finish
  %finish
%end

%external %real %fn %spec sin %alias "sinf" (%real angle)  ;! use C's math lib to check results.
%external %real %fn %spec cos %alias "cosf" (%real angle)

%real Theta
%integer i, sinTheta, cosTheta, CheckSinTheta, CheckCosTheta, sinErr = 0, cosErr = 0, err

%routine pround(%integer f, n)
  ! This is a quick (and rather ugly) hack due to the current Imp compiler
  ! not implementing all of the usual floating point I/O library calls yet
  
  %real r, sign = 1.0
  %integer tens = 0, p

  %if f=0 %start
    r = 0
  %else
    f = -f %and sign = -1.0 %if f < 0
    r = f / K16
    r = r * 10.0 %and tens = tens+1 %until 100 <= INT PT(r) < 1000
    p = tens
    r = INT(r)  { round to 3 significant digits }
    r = r / 10.0 %and tens = tens-1 %while tens > 0
    r = r + 0.000001 { handle real values not exactly representable in decimal }
  %finish
  spaces(n - (p+3))
  %if sign < 0 %then print symbol('-') %else print symbol(' ')
  print(r, p+2) { library call truncates, doesn't round :-( }
%end

  ! We get about 3 significant digits of accuracy
  print string("Note: this CORDIC was implemented to only 14 bits of accuracy for 16 bit machines,"); newline
  print string("      so we will report results to 3 significant digits."); newlines(2)
  
  SinCos(INT(-9*K64/(2*PI)), sinTheta, cosTheta)
  printstring("sin and cos of  -9 is "); pround(sinTheta,6); space; pround(cosTheta,6); newline

  SinCos(INT(0*K64/(2*PI)), sinTheta, cosTheta)
  printstring("sin and cos of   0 is "); pround(sinTheta,6); space; pround(cosTheta,6); newline

  SinCos(INT(1.5*K64/(2*PI)), sinTheta, cosTheta)
  printstring("sin and cos of 1.5 is "); pround(sinTheta,6); space; pround(cosTheta,6); newline

  SinCos(INT(6*K64/(2*PI)), sinTheta, cosTheta)
  printstring("sin and cos of   6 is "); pround(sinTheta,6); space; pround(cosTheta,6); newline


  %if debug # 0 %start
    %for i = 0, 4096, K64 %cycle  { A circle is 64K units }
      Theta = i*2*PI/K64 { A circle is 2Pi radians }
  
      SinCos(i, sinTheta, cosTheta)
  
      print string("Cordic:         Theta = "); phex(i)
      print string("; Sin(Theta) = "); write(sinTheta, 0)
      print string("; Cos(Theta) = "); write(cosTheta, 0)
      newline;
      
      CheckSinTheta = INT(sin(Theta)*K16); CheckCosTheta = INT(cos(Theta)*K16);
      
      print string("Floating point: Theta = "); phex(INT(Theta/(2*PI)*K64))
      print string("; Sin(Theta) = "); write(CheckSinTheta, 0)    
      print string("; Cos(Theta) = "); write(CheckCosTheta, 0)
      newlines(2)
  
      err = CheckSinTheta - sinTheta; err = -err %if err < 0; sinErr = sinErr + err
      err = CheckCosTheta - cosTheta; err = -err %if err < 0; cosErr = cosErr + err
  
    %repeat
  
    { Use this to test if algorithmic tweaks improve the accuracy }
    print string("Cumulative error:  sin = "); write(sinErr, 0)
    print string(",  cos = "); write(cosErr, 0)
    newline
    
  %finish
%end %of %program

Output:

Note: this CORDIC was implemented to only 14 bits of accuracy for 16 bit machines,
      so we will report results to 3 significant digits.

sin and cos of  -9 is -0.412 -0.911
sin and cos of   0 is  0.000   1.00
sin and cos of 1.5 is  0.998  0.0707
sin and cos of   6 is -0.280  0.960

J

Model implementation:

epsilon=: 10^-12
phin=: (#~ epsilon <: ])~.@(, _3 o. 10^-@#)^:_ ''
tent=: 10^-i.#phin

cordic=: {{alpha=. y
  XY=. 1 0 assert. 0 <: alpha assert. 0.25p1 >: alpha
  while. epsilon < alpha do.
    k=. phin I. alpha
    alpha=. alpha - k{phin
    XY =. XY + (k{tent) * XY +/ .* ((+.0j1 _1))
  end.
  XY
}}

CORDIC=: {{
  'octant angle'=. 8 0.25p1#:y
  select. octant
    case. 0 do.         cordic        angle
    case. 1 do.       |.cordic 0.25p1-angle
    case. 2 do. _1  1*|.cordic        angle
    case. 3 do. _1  1*  cordic 0.25p1-angle
    case. 4 do. _1 _1*  cordic        angle
    case. 5 do. _1 _1*|.cordic 0.25p1-angle
    case. 6 do.  1 _1*|.cordic        angle
    case. 7 do.  1 _1*  cordic 0.25p1-angle
  end.
}}

Task examples (cos is the left value in the result, argument is in radians):

   CORDIC -9
_0.92954 _0.420445
   CORDIC 0
1 0
   CORDIC 1.5
0.070762 0.997844
   CORDIC 6
0.970161 _0.282323

Notes

CAUTION: At the time of this writing, the task description declares that the cordic algorithm is valid in the range 0 .. π/2. But it appears that the algorithm can only be valid in the range 0..π/4.

In this J implementation, we use three constants, two of which are lookup tables (tent gives us negative powers of ten as a lookup table):

   epsilon
1e_12
   phin
0.785398 0.0996687 0.00999967 0.001 0.0001 1e_5 1e_6 1e_7 1e_8 1e_9 1e_10 1e_11 1e_12
   tent
1 0.1 0.01 0.001 0.0001 1e_5 1e_6 1e_7 1e_8 1e_9 1e_10 1e_11 1e_12

By default, J displays the first six digits of floating point numbers (floating point numbers are more precise, but in most circumstances the values being operated on are not more accurate than six digits). But of course, J retains the values at higher precision. For example:

   {&phin
{ &0.785398163397448279 0.0996686524911620381 0.00999966668666523936 0.000999999666666867007 9.99999996666667089e_5 9.99999999966667302e_6 9.99999999999667283e_7 9.9999999999999744e_8 1.00000000000000085e_8 1.00000000000000089e_9 1.00000000000000107e_10 1....

That said, it's not clear that this algorithm can be more accurate than something near 3.5% for the general case.

Also: double parenthesis around a noun phrase tells the interpreter that that expression is a constant which should be evaluated once, ahead of time.

Java

Translation of: C

public class CordicCalculator {

    // Constants
    private static final double[] angles = {
        0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
        0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
        0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119,
        0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812,
        0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863,
        0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929,
        0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058
    };

    private static final double[] kvalues = {
        0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775,
        0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889,
        0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894,
        0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314,
        0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925,
        0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888
    };

    // Convert degrees to radians
    private static double radians(double degrees) {
        return degrees * Math.PI / 180.0;
    }

    // Cordic algorithm implementation
    private static void cordic(double alpha, int n, double[] result) {
        int i, ix, sigma;
        double kn, x, y, atn, t, theta = 0.0, pow2 = 1.0;
        int newsgn = (int)Math.floor(alpha / (2.0 * Math.PI)) % 2 == 1 ? 1 : -1;
        if (alpha < -Math.PI/2.0 || alpha > Math.PI/2.0) {
            if (alpha < 0) {
                cordic(alpha + Math.PI, n, result);
            } else {
                cordic(alpha - Math.PI, n, result);
            }
            result[0] = result[0] * newsgn;
            result[1] = result[1] * newsgn;
            return;
        }
        ix = n - 1;
        if (ix > 23) ix = 23;
        kn = kvalues[ix];
        x = 1;
        y = 0;
        for (i = 0; i < n; ++i) {
            atn = angles[i];
            sigma = (theta < alpha) ? 1 : -1;
            theta += sigma * atn;
            t = x;
            x -= sigma * y * pow2;
            y += sigma * t * pow2;
            pow2 /= 2.0;
        }
        result[0] = x * kn;
        result[1] = y * kn;
    }

    // Main method
    public static void main(String[] args) {
        double c_cos, c_sin;
        double[] result = new double[2];
        double[] angles = {-9.0, 0.0, 1.5, 6.0};
        String format;
        
        System.out.println("  x       sin(x)     diff. sine     cos(x)    diff. cosine");
        format = "%+03d.0°  %+.8f (%+.8f) %+.8f (%+.8f)\n";
        for (int th = -90; th <= +90; th += 15) {
            double thr = radians(th);
            cordic(thr, 24, result);
            c_cos = result[0];
            c_sin = result[1];
            System.out.printf(format, th, c_sin, c_sin - Math.sin(thr), c_cos, c_cos - Math.cos(thr));
        }

        System.out.println("\nx(rads)   sin(x)     diff. sine     cos(x)    diff. cosine");
        format = "%+4.1f    %+.8f (%+.8f) %+.8f (%+.8f)\n";
        for (int i = 0; i < angles.length; ++i) {
            double thr = angles[i];
            cordic(thr, 24, result);
            c_cos = result[0];
            c_sin = result[1];
            System.out.printf(format, thr, c_sin, c_sin - Math.sin(thr), c_cos, c_cos - Math.cos(thr));
        }
    }
}

Output:

  x       sin(x)     diff. sine     cos(x)    diff. cosine
-90.0°  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
-75.0°  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
-60.0°  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
-45.0°  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
-30.0°  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
-15.0°  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
+00.0°  +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+15.0°  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
+30.0°  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
+45.0°  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
+60.0°  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
+75.0°  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
+90.0°  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

x(rads)   sin(x)     diff. sine     cos(x)    diff. cosine
-9.0    -0.41211842 (+0.00000006) -0.91113029 (-0.00000003)
+0.0    +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+1.5    +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
+6.0    -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)

jq

Adapted from Wren

Works with: jq

Works with gojq, the Go implementation of jq

Preliminaries

# like while/2 but emit the final term rather than the first one
def whilst(cond; update):
  def _whilst:
    if cond then update | (., _whilst) else empty end;
  _whilst;

# Simplistic approach to rounding:
def round($ndec): pow(10;$ndec) as $p | . * $p | round / $p;

def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;

The task

# The following are pre-computed to avoid using atan and sqrt functions.
def angles: [
    0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
    0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
    0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119,
    0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812,
    0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863,
    0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929,
    0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058
];

def kvalues: [
    0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775,
    0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889,
    0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894,
    0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314,
    0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925,
    0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888
];

def PI: (1 | atan * 4);

def radians: . * PI / 180;

def Cordic($alpha; $n):
  PI as $PI
  | (if (($alpha / (2 * $PI))|floor % 2 == 1) then 1 else -1 end) as $newsgn
  | if ($alpha < -$PI/2) 
    then Cordic($alpha + $PI; $n) | map(. * $newsgn)
    elif ($alpha >  $PI/2)
    then Cordic($alpha - $PI; $n) | map(. * $newsgn)
    else kvalues[ [$n-1, (kvalues|length-1)] | min] as $kn
    | reduce angles[0:$n][] as $atan ({ theta: 0, x: 1, y: 0, pow2: 1};
        (if .theta < $alpha then 1 else -1 end) as $sigma
        | .theta += $sigma * $atan
        | .x as $t
        | .x = (.x - $sigma * .y * .pow2)
        | .y = (.y + $sigma * $t * .pow2)
        | .pow2 /= 2 )
    | [.x * $kn, .y * $kn]
    end;

def example:
  def format: map(round(8) | lpad(11)) | join(" ");
  "    x         sin(x)     diff. sine   cos(x)   diff. cosine",
  ({ th: -90 }
   | whilst (.th <= 90;
       (.th | radians) as $thr
       | Cordic($thr; 24) as [ $cos, $sin ]
       | .out = [.th, $sin, $sin - ($thr|sin), $cos, $cos - ($thr|cos) ]
       | .th += 15 )
   | .out
   | format),

  "\n    x(rads)   sin(x)     diff. sine   cos(x)   diff. cosine",
  ((-9, 0, 1.5, 6) as $thr
   | Cordic($thr; 24) as [$cos, $sin]
   | [$thr, $sin, $sin - ($thr|sin), $cos, $cos - ($thr|cos)]
   | format) ;

example

Invocation: jq -ncr -f cordic.jq

Output:

    x         sin(x)     diff. sine   cos(x)   diff. cosine
        -90          -1           0      -7e-08      -7e-08
        -75 -0.96592585      -3e-08  0.25881895      -9e-08
        -60 -0.86602545      -5e-08  0.49999992      -8e-08
        -45 -0.70710684      -6e-08  0.70710672      -6e-08
        -30 -0.49999992       8e-08  0.86602545       5e-08
        -15 -0.25881895       9e-08  0.96592585       3e-08
          0       7e-08       7e-08           1          -0
         15  0.25881895      -9e-08  0.96592585       3e-08
         30  0.49999992      -8e-08  0.86602545       5e-08
         45  0.70710684       6e-08  0.70710672      -6e-08
         60  0.86602545       5e-08  0.49999992      -8e-08
         75  0.96592585       3e-08  0.25881895      -9e-08
         90           1          -0      -7e-08      -7e-08

    x(rads)   sin(x)     diff. sine   cos(x)   diff. cosine
         -9 -0.41211842       6e-08 -0.91113029      -3e-08
          0       7e-08       7e-08           1          -0
        1.5  0.99749499           0  0.07073719      -2e-08
          6 -0.27941552      -2e-08  0.96017028      -1e-08

Julia

""" Modified from MATLAB example code at en.wikipedia.org/wiki/CORDIC """

using Printf

"""
    Compute v = [cos(alpha), sin(alpha)] (alpha in radians).
    Increasing the iteration value will increase the precision.
"""
function cordic(alpha, iteration = 24)
    # Fix for the Wikipedia's MATLAB code bug in sin (sometimes flips sign) when |θ| > 2π
    newsgn = isodd(Int(floor(alpha / 2π))) ? 1 : -1
    alpha < -π/2 && return newsgn * cordic(alpha + π, iteration)
    alpha > π/2 && return newsgn * cordic(alpha - π, iteration)

    # Initialization of tables of constants used by CORDIC
    # need a table of arctangents of negative powers of two, in radians:
    # angles = atan(2.^-(0:27));
    angles =  [
        0.78539816339745,   0.46364760900081,   0.24497866312686,   0.12435499454676,
        0.06241880999596,   0.03123983343027,   0.01562372862048,   0.00781234106010,
        0.00390623013197,   0.00195312251648,   0.00097656218956,   0.00048828121119,
        0.00024414062015,   0.00012207031189,   0.00006103515617,   0.00003051757812,
        0.00001525878906,   0.00000762939453,   0.00000381469727,   0.00000190734863,
        0.00000095367432,   0.00000047683716,   0.00000023841858,   0.00000011920929,
        0.00000005960464,   0.00000002980232,   0.00000001490116,   0.00000000745058, ]
    # and a table of products of reciprocal lengths of vectors [1, 2^-2j]:
    # Kvalues = cumprod(1./sqrt(1 + 1j*2.^(-(0:23))))
    Kvalues = [
        0.70710678118655,   0.63245553203368,   0.61357199107790,   0.60883391251775,
        0.60764825625617,   0.60735177014130,   0.60727764409353,   0.60725911229889,
        0.60725447933256,   0.60725332108988,   0.60725303152913,   0.60725295913894,
        0.60725294104140,   0.60725293651701,   0.60725293538591,   0.60725293510314,
        0.60725293503245,   0.60725293501477,   0.60725293501035,   0.60725293500925,
        0.60725293500897,   0.60725293500890,   0.60725293500889,   0.60725293500888, ]

    Kn = Kvalues[min(iteration, length(Kvalues))]  
    # Initialize loop variables:
    v = [1, 0] # start with 2-vector cosine and sine of zero
    poweroftwo = 1
    angle = angles[1]
    
    # Iterations
    for j = 0:iteration-1
        if alpha < 0
            sigma = -1
        else
            sigma = 1
        end
        factor = sigma * poweroftwo
        # Note the matrix multiplication can be done using scaling by powers of two and addition subtraction
        R = [1 -factor
             factor  1]
        v = R * v # 2-by-2 matrix multiply
        alpha -= sigma * angle # update the remaining angle
        poweroftwo /= 2
        # update the angle from table, or eventually by just dividing by two
        if j + 2 > length(angles)
            angle /= 2
        else
            angle = angles[j + 2]
        end
    end
    
    # Adjust length of output vector to be [cos(alpha), sin(alpha)]:
    v .*= Kn
    return v
end

function test_cordic()
    println("  x       sin(x)     diff. sine     cos(x)    diff. cosine ")
    for θ in -90:15:90 
        cosθ, sinθ = cordic(deg2rad(θ))
        @printf("%+05.1f°  %+.8f (%+.8f) %+.8f (%+.8f)\n",
           θ, sinθ, sinθ - sind(θ), cosθ, cosθ - cosd(θ))
    end
    println("\nx(radians)  sin(x)     diff. sine     cos(x)    diff. cosine ")
    for θr in [-9, 0, 1.5, 6]
        cosθ, sinθ = cordic(θr)
        @printf("%+3.1f      %+.8f (%+.8f) %+.8f (%+.8f)\n",
           θr, sinθ, sinθ - sin(θr), cosθ, cosθ - cos(θr))
    end

end

test_cordic()

Output:

  x       sin(x)     diff. sine     cos(x)    diff. cosine 
-90.0°  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
-75.0°  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
-60.0°  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
-45.0°  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
-30.0°  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
-15.0°  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
+00.0°  -0.00000007 (-0.00000007) +1.00000000 (-0.00000000)
+15.0°  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
+30.0°  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
+45.0°  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
+60.0°  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
+75.0°  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
+90.0°  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

x(radians)  sin(x)     diff. sine     cos(x)    diff. cosine
-9.0      -0.41211842 (+0.00000006) -0.91113029 (-0.00000003)
+0.0      -0.00000007 (-0.00000007) +1.00000000 (-0.00000000)
+1.5      +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
+6.0      -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)

Lua

Translation of: ALGOL 68

do -- implement sin, cos and tan using the CORDIC algorithm

    local piBy2 = math.pi / 2

    -- pre-computed table of arctan(10^-n) values
    local theta = { 7.85398163397448e-01, 9.96686524911620e-02
                  , 9.99966668666524e-03, 9.99999666666867e-04
                  , 9.99999996666667e-05, 9.99999999966667e-06
                  , 9.99999999999667e-07, 9.99999999999997e-08
                  , 1.00000000000000e-08, 1.00000000000000e-09
                  , 1.00000000000000e-10, 1.00000000000000e-11
                  , 1.00000000000000e-12, 1.00000000000000e-13
                  , 1.00000000000000e-14, 1.00000000000000e-15
                  , 1.00000000000000e-16
                  }
    local epsilon = 1e-16

    --[[ CORDIC algorithm, finds "y" and "x" for alpha radians
         signs indicates the sign of the results:
               signs[ 1 ] = sign for -π : -π/2, signs[ 1 ] = sign for -π/2 : 0
               signs[ 3 ] = sign for  0 :  π/2, signs[ 4 ] = sign for  π/2 : π
         signBoth = TRUE => sign applied to both y and x, FALSE => sign y only
    --]]
    function cordic( alphaIn, signs, signBoth )
        local alpha     = alphaIn
        -- ensure -π <= alpha <= π
        local flipSign = false
        while alpha < - math.pi do alpha = alpha + math.pi flipSign = not flipSign end
        while alpha >   math.pi do alpha = alpha - math.pi flipSign = not flipSign end
        local sign
        if    alpha < - piBy2 then
            alpha = alpha + math.pi
            sign = signs[ 1 ]
        elseif  alpha < 0     then
            alpha = - alpha
            sign = signs[ 2 ]
        elseif  alpha < piBy2 then
            sign = signs[ 3 ]
        else -- alpha <= math.pi
            alpha = math.pi - alpha
            sign = signs[ 4 ]
        end
        if flipSign and signBoth then sign = -sign end
        local x, y, k, tenToMinusK = 1, 0, 1, 1 -- NB: tenToMinusK is 10^-(k-1)
        while epsilon < alpha do
            while alpha < theta[ k ] do
                k           = k +  1
                tenToMinusK = tenToMinusK / 10
            end
            alpha = alpha - theta[ k ]
            x, y = x - tenToMinusK * y, y + tenToMinusK * x
        end
        return sign * y, signBoth and sign * x or x
    end

    -- cos(alpha) using the CORDIC algorithm. alpha in radians
    function cCos( alpha )
        local y, x = cordic( alpha, { -1, 1, 1, -1 }, true )
        return x / math.sqrt( ( x * x ) + ( y * y ) )
    end

    -- sin(alpha) using the CORDIC algorithm. alpha in radians
    function cSin( alpha )
        local y, x = cordic( alpha, { -1,-1, 1, 1 }, true )
        return y / math.sqrt( ( x * x ) + ( y * y ) )
    end

    -- tan(alpha) using the CORDIC algorithm. alpha in radians
    function cTan( alpha )
        local y, x  = cordic( alpha, {  1, -1, 1,-1 }, false )
        return x == 0 and NaN or y / x
    end

    function f6( v ) return string.format( " %9.6f", v ) end
    function g(  v ) return string.format( " %12g",  v ) end

    function showCordic( angle )
        local cosine,  cordicCos = math.cos( angle ), cCos( angle )
        local sine,    cordicSin = math.sin( angle ), cSin( angle )
        local tangent, cordicTan = math.tan( angle ), cTan( angle )
        local cDiff  = cordicCos - cosine  if cDiff < 0 then cDiff = - cDiff end
        local sDiff  = cordicSin - sine    if sDiff < 0 then sDiff = - sDiff end
        local tDiff  = cordicTan - tangent if tDiff < 0 then tDiff = - tDiff end
        io.write( string.format( "%4.1f: ", angle )
                , f6( cordicCos ), f6( cosine  ), g( cDiff ), " |"
                , f6( cordicSin ), f6( sine    ), g( sDiff ), " |"
                , f6( cordicTan ), f6( tangent ), g( tDiff ), "\n"
                )
    end 

    local tests = { -9, 0, 1.5, 6 }
    io.write( "angle cordic cos       cos   difference" )
    io.write(     "  cordic sin       sin   difference" )
    io.write(     "  cordic tan       tan   difference" )
    io.write( "\n" )
    for i = 1, # tests do
        showCordic( tests[ i ] )
    end
    io.write( "angle cordic cos       cos   difference" )
    io.write(     "  cordic sin       sin   difference" )
    io.write(     "  cordic tan       tan   difference" )
    io.write( "\n" )
    for i = 7810, 7891, 9 do
        showCordic( i / 10000 )
    end

end

Output:

angle cordic cos       cos   difference  cordic sin       sin   difference  cordic tan       tan   difference
-9.0:  -0.911130 -0.911130  1.11022e-16 | -0.412118 -0.412118            0 |  0.452316  0.452316            0
 0.0:   1.000000  1.000000            0 |  0.000000  0.000000            0 |  0.000000  0.000000            0
 1.5:   0.070737  0.070737  7.49401e-16 |  0.997495  0.997495            0 | 14.101420 14.101420   1.5099e-13
 6.0:   0.960170  0.960170  1.11022e-16 | -0.279415 -0.279415  5.55112e-16 | -0.291006 -0.291006  6.10623e-16
angle cordic cos       cos   difference  cordic sin       sin   difference  cordic tan       tan   difference
 0.8:   0.710210  0.710210  3.33067e-16 |  0.703990  0.703990  2.22045e-16 |  0.991242  0.991242  6.66134e-16
 0.8:   0.709576  0.709576  2.22045e-16 |  0.704629  0.704629  1.11022e-16 |  0.993028  0.993028  4.44089e-16
 0.8:   0.708942  0.708942  3.33067e-16 |  0.705267  0.705267  4.44089e-16 |  0.994817  0.994817  9.99201e-16
 0.8:   0.708307  0.708307  6.66134e-16 |  0.705905  0.705905  5.55112e-16 |  0.996609  0.996609  1.66533e-15
 0.8:   0.707671  0.707671  5.55112e-16 |  0.706542  0.706542  3.33067e-16 |  0.998405  0.998405  1.33227e-15
 0.8:   0.707035  0.707035  1.11022e-16 |  0.707179  0.707179  1.11022e-16 |  1.000204  1.000204  2.22045e-16
 0.8:   0.706398  0.706398  1.11022e-16 |  0.707815  0.707815  2.22045e-16 |  1.002006  1.002006  4.44089e-16
 0.8:   0.705761  0.705761  1.11022e-16 |  0.708450  0.708450  1.11022e-16 |  1.003811  1.003811  2.22045e-16
 0.8:   0.705123  0.705123  3.33067e-16 |  0.709085  0.709085  5.55112e-16 |  1.005619  1.005619  1.33227e-15
 0.8:   0.704484  0.704484  5.55112e-16 |  0.709720  0.709720  4.44089e-16 |  1.007431  1.007431  1.33227e-15

Perl

Translation of: Raku

use strict;
use warnings;

sub CORDIC {
    my($a) = shift;
    my ($k, $x, $y) = (0, 1, 0);
    my @PoT = (1, 0.1, 0.01, 0.001, 0.0001, 0.00001);
    my @Tbl = <7.853981633e-1 9.966865249e-2 9.999666686e-3 9.999996666e-4 9.999999966e-5 9.999999999e-6 0>;

    while ($a > 1e-5) {
        $k++ while $a < $Tbl[$k];
        $a -= $Tbl[$k];
        ($x,$y) = ($x - $PoT[$k]*$y, $y + $PoT[$k]*$x);
    }
    $x / sqrt($x*$x + $y*$y)
}

print "Angle    CORDIC       Cosine       Error\n";
for my $angle (<-9 0 1.5 6>) {
    my $cordic = CORDIC abs $angle;
    printf "%4.1f %12.8f %12.8f %12.8f\n", $angle, $cordic, cos($angle), cos($angle) - $cordic
}

Output:

Angle    CORDIC       Cosine       Error
-9.0  -0.91112769  -0.91113026  -0.00000257
 0.0   1.00000000   1.00000000   0.00000000
 1.5   0.07073880   0.07073720  -0.00000160
 6.0   0.96016761   0.96017029   0.00000268

Phix

Translation of: Wren

with javascript_semantics
constant angles = {
    0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
    0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
    0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119,
    0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812,
    0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863,
    0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929,
    0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058}
constant kvalues = {
    0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775,
    0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889,
    0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894,
    0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314,
    0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925,
    0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888}

function cordic(atom alpha, integer n)
    bool sgn = +1
    while alpha < -PI/2 do alpha += PI sgn *= -1 end while
    while alpha >  PI/2 do alpha -= PI sgn *= -1 end while
    atom kn = sgn * kvalues[min(n,length(kvalues))],
         atan, theta = 0, x = 1, y = 0, pow2 = 1
    for i=1 to n do
        atan = iff(i<=length(angles)?angles[i]:atan/2)
        atom sigma = iff(theta < alpha ? 1 : -1), t = x
--      atom sigma = iff(theta <= alpha ? 1 : -1), t = x -- (matches Julia)
        theta += sigma * atan
        x -= sigma * y * pow2
        y += sigma * t * pow2
        pow2 /= 2
    end for
    return {x * kn, y * kn}
end function

constant fmh = "%n x(%s)  sin(x)       diff.sine    cos(x)      diff.cosine\n",
         fmt = "%+6.1f  %+.8f (%+.8f) %+.8f (%+.8f)\n"
printf(1,fmh,{false,"deg"})
for theta = -90 to 90 by 15 do
    atom r = theta*PI/180,
       {c,s} = cordic(r,24)
    printf(1,fmt,{theta,s,s-sin(r),c,c-cos(r)})
end for
printf(1,fmh,{true,"rad"})
for r in {-9, 0, 1.5, 6} do
    atom {c,s} = cordic(r, 24)
    printf(1,fmt,{r,s,s-sin(r),c,c-cos(r)})
end for

Output:

As noted above in the commented out line, and mentioned in the Wren entry, testing for <= flips the sign of sin(0).

 x(deg)  sin(x)       diff.sine    cos(x)      diff.cosine
 -90.0  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
 -75.0  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
 -60.0  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
 -45.0  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
 -30.0  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
 -15.0  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
  +0.0  +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
 +15.0  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
 +30.0  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
 +45.0  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
 +60.0  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
 +75.0  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
 +90.0  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

 x(rad)  sin(x)       diff.sine    cos(x)      diff.cosine
  -9.0  -0.41211842 (+0.00000006) -0.91113029 (-0.00000003)
  +0.0  +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
  +1.5  +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
  +6.0  -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)

Python

Translation of: C

# CORDIC.py by Xing216
from math import pi, sin, cos, floor

angles = [
    0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
    0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
    0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119,
    0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812,
    0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863,
    0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929,
    0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058
]
kvalues = [
    0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775,
    0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889,
    0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894,
    0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314,
    0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925,
    0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888
]
radians = lambda degrees: (degrees * pi) / 180

def Cordic(alpha, n, c_cos, c_sin):
    i, ix, sigma = 0, 0, 0
    kn, x, y, atn, t, theta = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
    pow2 = 1.0
    newsgn = 1 if floor(alpha / (2.0 * pi)) % 2 == 1 else -1
    if alpha < -pi/2.0 or alpha > pi/2.0:
        if alpha < 0:
            Cordic(alpha + pi, n, x, y)
        else:
            Cordic(alpha - pi, n, x, y)
        c_cos = x * newsgn
        c_sin = y * newsgn
        return c_cos, c_sin
    ix = n - 1
    if ix > 23: ix = 23
    kn = kvalues[ix]
    x = 1
    y = 0
    for i in range(n):
        atn = angles[i]
        sigma = 1 if theta < alpha else -1
        theta += sigma * atn
        t = x
        x -= sigma * y * pow2
        y += sigma * t * pow2
        pow2 /= 2.0
    c_cos = x * kn
    c_sin = y * kn
    return c_cos, c_sin

def main():
    i, th = 0, 0
    thr, c_cos, c_sin = 0.0, 0.0, 0.0
    angles = [-9.0, 0.0, 1.5, 6.0]
    print("  x       sin(x)     diff. sine     cos(x)    diff. cosine")
    for th in range(-90,91,15):
        thr = radians(th)
        c_cos, c_sin = Cordic(thr, 24, c_cos, c_sin)
        sin_diff = c_sin - sin(thr)
        cos_diff = c_cos - cos(thr)
        print(f"{th:+03}.0°  {c_sin:+.8f} ({sin_diff:+.8f}) {c_cos:+.8f} ({cos_diff:+.8f})")
    print("\nx(rads)   sin(x)     diff. sine     cos(x)    diff. cosine")
    for i in range(4):
        thr = angles[i]
        c_cos, c_sin = Cordic(thr, 24, c_cos, c_sin)
        sin_diff = c_sin - sin(thr)
        cos_diff = c_cos - cos(thr)
        print(f"{thr:+4.1f}    {c_sin:+.8f} ({c_sin - sin(thr):+.8f}) {c_cos:+.8f} ({c_cos - cos(thr):+.8f})")

if __name__ in "__main__":
    main()

Output:

  x       sin(x)     diff. sine     cos(x)    diff. cosine
-90.0°  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
-75.0°  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
-60.0°  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
-45.0°  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
-30.0°  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
-15.0°  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
+00.0°  +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+15.0°  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
+30.0°  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
+45.0°  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
+60.0°  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
+75.0°  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
+90.0°  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

x(rads)   sin(x)     diff. sine     cos(x)    diff. cosine
-9.0    -0.00000000 (+0.41211849) -0.00000000 (+0.91113026)
+0.0    +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+1.5    +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
+6.0    -0.00000000 (+0.27941550) -0.00000000 (-0.96017029)

Raku

Translation of: XPL0

# 2023082 Raku programming solution

sub CORDIC ($A is copy) {

   my (\Ten, $K, $X, $Y) = ( 1, * * 1/10 ... * )[^6], 0, 1, 0;

   my \Tbl = < 7.853981633974480e-1 9.966865249116200e-2 9.999666686665240e-3 
               9.999996666668670e-4 9.999999966666670e-5 9.999999999666670e-6 0.0>;

   while $A > 1e-5 { 
      $K++ while $A < Tbl[$K];
      $A -= Tbl[$K];
      ($X,$Y) = $X - Ten[$K]*$Y, $Y + Ten[$K]*$X;
   }
   return $X, sqrt($X*$X + $Y*$Y)
}

say "Angle    CORDIC       Cosine       Error";
for <-9 0 1.5 6> { 
   my \result = [/] CORDIC .abs;
   printf "% 2.1f  "~"% 2.8f  " x 3~"\n", $_, result, .cos, .cos - result
}

Output:

Same as XPL0 or you may Attempt This Online!

RPL

Works with: HP version 28

≪ RAD { } 1
   DO
      SWAP OVER ATAN + SWAP 10 /
   UNTIL DUP DUP TAN == END        @ memorize constants until precision limit is reached 
   DROP 'THN' STO 
   THN SIZE →STR " constants in memory." *
   1E-12 'EPSILON' STO  
≫ ≫ 'INIT' STO       

≪ IF DUP THEN 
      1 SWAP START 10 / NEXT       @ shift one digit right 
   ELSE DROP END
≫ 'SR10' STO 

≪ IF THN SIZE OVER 1 + <    
   THEN 1 SWAP SR10                @ get arctan(θ[k]) from memory 
   ELSE THN SWAP 1 + GET END       @ arctan(θ[k]) ≈ θ[k]
≫ '→THK' STO 

≪ → alpha
  ≪ 0 1 0                         @ initialize y, x and k
     WHILE alpha EPSILON > REPEAT
        WHILE DUP →THK alpha > REPEAT 
           1 + END
        'alpha' OVER →THK STO-
        DUP2 SR10 4 PICK + 4 ROLLD
        ROT OVER SR10 ROT SWAP -
        SWAP
     END 
     DROP /
≫ ≫ '→TAN' STO    

≪ 1 CF 
   '2*π' →NUM MOD
   IF DUP π / 2 * →NUM IP THEN
      { ≪ π SWAP 1 SF ≫ ≪ π 1 SF ≫ ≪ '2*π' SWAP ≫ }    @ corrections for angles > π/2
      LASTARG GET EVAL →NUM -                              @ apply correction according to quadrant
      END
   →TAN SQ 1 + √ INV
   IF 1 FS? THEN NEG END
≫ '→COS' STO  
     
≪ INIT { -9 0 1.5 6 } { }
   1 3 PICK SIZE FOR j 
      OVER j GET →COS + 
   NEXT SWAP DROP
≫ 'TASK' STO

Output:

2: "6 constants in memory."
1: { -.91113026188 1 7.07372016661E-2 .960170286655 }

Rust

Translation of: C

use std::f64::consts::PI;

// Constants for angles and kvalues
const ANGLES: [f64; 28] = [
    0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
    0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
    0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119,
    0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812,
    0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863,
    0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929,
    0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058,
];

const KVALUES: [f64; 24] = [
    0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775,
    0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889,
    0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894,
    0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314,
    0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925,
    0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888,
];

// Function to convert degrees to radians
fn radians(degrees: f64) -> f64 {
    degrees * PI / 180.0
}

// Cordic algorithm implementation
fn cordic(alpha: f64, n: usize, c_cos: &mut f64, c_sin: &mut f64) {
    let mut theta = 0.0;
    let mut pow2 = 1.0;
    let mut x = 1.0;
    let mut y = 0.0;
    let newsgn = if (alpha / (2.0 * PI)).floor() as i32 % 2 == 1 { 1.0 } else { -1.0 };

    if alpha < -PI / 2.0 || alpha > PI / 2.0 {
        if alpha < 0.0 {
            cordic(alpha + PI, n, &mut x, &mut y);
        } else {
            cordic(alpha - PI, n, &mut x, &mut y);
        }
        *c_cos = x * newsgn;
        *c_sin = y * newsgn;
        return;
    }

    let ix = if n - 1 > 23 { 23 } else { n - 1 };
    let kn = KVALUES[ix];

    for i in 0..n {
        let atn = ANGLES[i];
        let sigma = if theta < alpha { 1.0 } else { -1.0 };
        theta += sigma * atn;
        let t = x;
        x -= sigma * y * pow2;
        y += sigma * t * pow2;
        pow2 /= 2.0;
    }

    *c_cos = x * kn;
    *c_sin = y * kn;
}

fn main() {
    let mut c_cos: f64;
    let mut c_sin: f64;
    let test_angles = [-9.0, 0.0, 1.5, 6.0];

    println!("  x       sin(x)     diff. sine     cos(x)    diff. cosine");
    for th in (-90..=90).step_by(15) {
        let thr = radians(th as f64);
        c_cos = 0.0;
        c_sin = 0.0;
        cordic(thr, 24, &mut c_cos, &mut c_sin);
        println!("{:+03}.0°  {:+.8} ({:+.8}) {:+.8} ({:+.8})", th, c_sin, c_sin - thr.sin(), c_cos, c_cos - thr.cos());
    }

    println!("\nx(rads)   sin(x)     diff. sine     cos(x)    diff. cosine");
    for &angle in &test_angles {
        let thr = angle;
        c_cos = 0.0;
        c_sin = 0.0;
        cordic(thr, 24, &mut c_cos, &mut c_sin);
        println!("{:+4.1}    {:+.8} ({:+.8}) {:+.8} ({:+.8})", thr, c_sin, c_sin - thr.sin(), c_cos, c_cos - thr.cos());
    }
}

Output:

  x       sin(x)     diff. sine     cos(x)    diff. cosine
-90.0°  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
-75.0°  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
-60.0°  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
-45.0°  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
-30.0°  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
-15.0°  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
+00.0°  +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+15.0°  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
+30.0°  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
+45.0°  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
+60.0°  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
+75.0°  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
+90.0°  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

x(rads)   sin(x)     diff. sine     cos(x)    diff. cosine
-9.0    -0.41211842 (+0.00000006) -0.91113029 (-0.00000003)
+0.0    +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+1.5    +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
+6.0    -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)

Sidef

Translation of: Perl

func CORDIC(a) {
    var (k, x, y) = (0, 1, 0)
    var PoT = %n(1, 0.1, 0.01, 0.001, 0.0001, 0.00001)
    var Tbl = %n(7.853981633e-1 9.966865249e-2 9.999666686e-3 9.999996666e-4 9.999999966e-5 9.999999999e-6 0)

    for (true; a > 1e-5; a -= Tbl[k]) {
        while (a < Tbl[k]) { ++k }
        (x,y) = (x - PoT[k]*y, y + PoT[k]*x)
    }
    x / hypot(x,y)
}

print("Angle    CORDIC       Cosine       Error\n")
for angle in (%n(-9 0 1.5 6)) {
    var cordic = CORDIC(angle.abs)
    printf("%4.1f %12.8f %12.8f %12.8f\n", angle, cordic, cos(angle), cos(angle) - cordic)
}

Output:

Angle    CORDIC       Cosine       Error
-9.0  -0.91112769  -0.91113026  -0.00000257
 0.0   1.00000000   1.00000000   0.00000000
 1.5   0.07073880   0.07073720  -0.00000160
 6.0   0.96016761   0.96017029   0.00000268

Wren

Library: Wren-fmt

This is based on the Python example in the Wikipedia article except that the constants have been pre-calculated and the angles are adjusted, where necessary, to bring them within the [-π/2, π/2] range. The number of iterations is taken as 24 to try and match the Julia output which it does except for sin(0) which, curiously, has the same magnitude but a different sign.

import "./fmt" for Fmt

// The following are pre-computed to avoid using atan and sqrt functions.
var angles = [
    0.78539816339745, 0.46364760900081, 0.24497866312686, 0.12435499454676,
    0.06241880999596, 0.03123983343027, 0.01562372862048, 0.00781234106010,
    0.00390623013197, 0.00195312251648, 0.00097656218956, 0.00048828121119,
    0.00024414062015, 0.00012207031189, 0.00006103515617, 0.00003051757812,
    0.00001525878906, 0.00000762939453, 0.00000381469727, 0.00000190734863,
    0.00000095367432, 0.00000047683716, 0.00000023841858, 0.00000011920929,
    0.00000005960464, 0.00000002980232, 0.00000001490116, 0.00000000745058
]

var kvalues = [
    0.70710678118655, 0.63245553203368, 0.61357199107790, 0.60883391251775,
    0.60764825625617, 0.60735177014130, 0.60727764409353, 0.60725911229889,
    0.60725447933256, 0.60725332108988, 0.60725303152913, 0.60725295913894,
    0.60725294104140, 0.60725293651701, 0.60725293538591, 0.60725293510314,
    0.60725293503245, 0.60725293501477, 0.60725293501035, 0.60725293500925,
    0.60725293500897, 0.60725293500890, 0.60725293500889, 0.60725293500888
]

var PI = Num.pi

var radians = Fn.new { |d| d * PI / 180 }

var Cordic = Fn.new { |alpha, n|
    var newsgn = ((alpha / (2 * PI)).floor % 2 == 1) ? 1 : -1
    if (alpha < -PI/2) {
        var res = Cordic.call(alpha + PI, n)
        return [newsgn * res[0], newsgn * res[1]]
    }
    if (alpha >  PI/2) {
        var res = Cordic.call(alpha - PI, n)
        return [newsgn * res[0], newsgn * res[1]]
    }
    var kn = kvalues[(n-1).min(kvalues.count-1)]
    var theta = 0
    var x = 1
    var y = 0
    var pow2 = 1
    for (atan in angles[0...n]) {
        var sigma = (theta < alpha) ? 1 : -1
        theta = theta + sigma * atan
        var t = x
        x = x - sigma * y * pow2
        y = y + sigma * t * pow2
        pow2 = pow2 / 2
    }
    return [x * kn, y * kn]
}

Fmt.print("  x       sin(x)     diff. sine     cos(x)    diff. cosine")
var f = "$+5.1f°  $+.8f ($+.8f) $+.8f ($+.8f)"

var th = -90
while (th <= 90) {
    var thr = radians.call(th)
    var res = Cordic.call(thr, 24)
    var cos = res[0]
    var sin = res[1]
    Fmt.print(f, th, sin, sin - thr.sin, cos, cos - thr.cos)
    th = th + 15
}
f = "$+5.1f°  $+.8f ($+.8f) $+.8f ($+.8f)"
Fmt.print("\nx(rads)   sin(x)     diff. sine     cos(x)    diff. cosine ")
f = "$+4.1f    $+.8f ($+.8f) $+.8f ($+.8f)"
for (thr in [-9, 0, 1.5, 6]) {
    var res = Cordic.call(thr, 24)
    var cos = res[0]
    var sin = res[1]
    Fmt.print(f, thr, sin, sin - thr.sin, cos, cos - thr.cos)
}

Output:

  x       sin(x)     diff. sine     cos(x)    diff. cosine
-90.0°  -1.00000000 (+0.00000000) -0.00000007 (-0.00000007)
-75.0°  -0.96592585 (-0.00000003) +0.25881895 (-0.00000009)
-60.0°  -0.86602545 (-0.00000005) +0.49999992 (-0.00000008)
-45.0°  -0.70710684 (-0.00000006) +0.70710672 (-0.00000006)
-30.0°  -0.49999992 (+0.00000008) +0.86602545 (+0.00000005)
-15.0°  -0.25881895 (+0.00000009) +0.96592585 (+0.00000003)
 +0.0°  +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+15.0°  +0.25881895 (-0.00000009) +0.96592585 (+0.00000003)
+30.0°  +0.49999992 (-0.00000008) +0.86602545 (+0.00000005)
+45.0°  +0.70710684 (+0.00000006) +0.70710672 (-0.00000006)
+60.0°  +0.86602545 (+0.00000005) +0.49999992 (-0.00000008)
+75.0°  +0.96592585 (+0.00000003) +0.25881895 (-0.00000009)
+90.0°  +1.00000000 (-0.00000000) -0.00000007 (-0.00000007)

x(rads)   sin(x)     diff. sine     cos(x)    diff. cosine 
-9.0    -0.41211842 (+0.00000006) -0.91113029 (-0.00000003)
+0.0    +0.00000007 (+0.00000007) +1.00000000 (-0.00000000)
+1.5    +0.99749499 (+0.00000000) +0.07073719 (-0.00000002)
+6.0    -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)

XPL0

include xpllib;         \for Print

real X, Y, R;

proc CORDIC(A);
real A;
real Tbl, Ten, T;
int  K;
[Ten:= [1E0, 1E-1, 1E-2, 1E-3, 1E-4, 1E-5];
Tbl:= [7.853981633974480E-001,
       9.966865249116200E-002,
       9.999666686665240E-003,
       9.999996666668670E-004,
       9.999999966666670E-005,
       9.999999999666670E-006,
       0.0];
X:= 1.;  Y:= 0.;  K:= 0;
while A > 1E-5 do
        [while A < Tbl(K) do K:= K+1;
        A:= A - Tbl(K);
        T:= X - Ten(K)*Y;
        Y:= Y + Ten(K)*X;
        X:= T;
        ];
R:= sqrt(X*X + Y*Y);
];

real Angles, A;
int  I;
[Print("Angle    CORDIC       Cosine       Error\n");
Angles:= [-9., 0., 1.5, 6.];
for I:= 0 to 3 do
        [A:= Angles(I);
        CORDIC(abs(A));
        Print("%2.1f  %2.8f  %2.8f  %2.8f\n", A, X/R, Cos(A), Cos(A)-X/R);
        ];
]

Output:

Angle    CORDIC       Cosine       Error
-9.0  -0.91112769  -0.91113026  -0.00000257
 0.0   1.00000000   1.00000000   0.00000000
 1.5   0.07073880   0.07073720  -0.00000160
 6.0   0.96016761   0.96017029   0.00000268