Rosetta Code

CORDIC

Revision as of 03:50, 9 July 2023 by Wherrera (talk | contribs) (→‎{{header|Julia}}: cleanup)

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.

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

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 precision < 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

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 cosine 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)

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 }
Retrieved from "https://rosettacode.org/wiki/CORDIC?oldid=346073"