CORDIC: Difference between revisions

From Rosetta Code
Content added Content deleted
Line 117: Line 117:


test_cordic()
test_cordic()
</syntaxhigh;ight>{{out}}
</syntaxhighlight>{{out}}
<pre>
<pre>
x sin(x) diff. sine cos(x) diff. cosine
x sin(x) diff. sine cos(x) diff. cosine

Revision as of 03:48, 9 July 2023

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 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 = vec([
        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=346072"