CORDIC
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
FreeBASIC
#define min(a, b) iif((a) < (b), (a), (b))
#define floor(x) ((x*2.0-0.5) Shr 1)
#define pi 4 * Atn(1)
#define radians(x) ((x) * pi / 180)
Function CORDIC(alfa As Integer, iteracion As Integer = 24) As Double
Dim As Double v
' This function computes v = [cos(alpha), sin(alpha)] (alpha in radians)
' using iteration increasing iteration value will increase the precision.
If alfa < -pi/2 Or alfa > pi/2 Then
v = Iif(alfa < 0, CORDIC(alfa + pi, iteracion), CORDIC(alfa - pi, iteracion))
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 angulos(1 To 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}
' 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 Kvalores(1 To 28) = { _
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 = Kvalores(min(iteracion, Len(Kvalores)))
' Initialize loop variables:
Dim As Integer poderde2 = 1, sigma, factor, R
Dim As Double angulo = angulos(1)
' iteracions
For j As Integer = 0 To iteracion-1
sigma = Iif(alfa < 0, -1, 1)
factor = sigma * poderde2
' Note the matrix multiplication can be done using scaling by powers of two and addition subtraction
R = factor
v *= R ' 2-by-2 matrix multiply
alfa -= sigma * angulo ' update the remaining angulo
poderde2 /= 2
' update the angulo from table, or eventually by just dividing by two
If j + 2 > Len(angulos) Then
angulo /= 2
Else
angulo = angulos(j + 2)
End If
Next
' Adjust length of output vector to be (cos(alfa), sin(alfa)):
v *= Kn
Return v
End Function
Dim As Double test(1 To 4) = {-9, 0, 1.5, 6}
Print !"\nx(radians) cos(x)"
For r As Integer = 1 To 4
Print Using " +#.# +#.########"; test(r); Cos(test(r))
Next
Sleep- Output:
x(radians) cos(x) -9.0 -0.91113026 +0.0 +1.00000000 +1.5 +0.07073720 +6.0 +0.96017029
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 2% 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.
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
≪ 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 }
Wren
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 res = Cordic.call(radians.call(th), 24)
var cos = res[0]
var sin = res[1]
var thr = radians.call(th)
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