CORDIC: Difference between revisions
m
→{{header|Phix}}: use pygments
m (→{{header|Julia}}: cleanup) |
m (→{{header|Phix}}: use pygments) |
||
| (39 intermediate revisions by 14 users not shown) | |||
Line 1:
{{
;Introduction
Line 15:
function tan(alpha) <span style="color:grey">// 0 < alpha <= π/2 </span>
x = 1 ; y = 0 ; k = 0
while
while alpha < θ[k]
k++
Line 30:
*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
<br>
;See also
[[wp:CORDIC|CORDIC on Wikipedia]]
<br>
=={{header|ALGOL 68}}==
Based on the pseudo code, method to calculate cos (and sin) from the XPL0 sample.
<syntaxhighlight lang="algol68">
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
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|C}}==
{{trans|Wren}}
If the M_PI constant is not defined in your ''math.h'', then you can use 4 * atan(1) instead.
<syntaxhighlight lang="c">#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;
}</syntaxhighlight>
{{out}}
<pre>
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)
</pre>
=={{header|FreeBASIC}}==
{{incorrect|FreeBASIC| <br><br> The provided function is not used. The answers provided are rather via call to a FreeBASIC library function.}}
{{trans|Julia}}
<syntaxhighlight lang="vb">#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</syntaxhighlight>
{{out}}
<pre>x(radians) cos(x)
-9.0 -0.91113026
+0.0 +1.00000000
+1.5 +0.07073720
+6.0 +0.96017029</pre>
=={{header|J}}==
Model implementation:
<syntaxhighlight lang=J>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.
}}</syntaxhighlight>
Task examples (cos is the left value in the result, argument is in radians):
<syntaxhighlight lang=J> CORDIC -9
_0.92954 _0.420445
CORDIC 0
1 0
CORDIC 1.5
0.070762 0.997844
CORDIC 6
0.970161 _0.282323</syntaxhighlight>
=== 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 (<code>tent</code> gives us negative powers of ten as a lookup table):
<syntaxhighlight lang=J> 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</syntaxhighlight>
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:
<syntaxhighlight lang=J> {&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....</syntaxhighlight>
That said, it's not clear that this algorithm can be more accurate than something near 3.5% for the general case.
Also: [[j:Vocabulary/SpecialCombinations#Preexecuting_verbs_with_.28.28_.29.29|double parenthesis around a noun phrase]] tells the interpreter that that expression is a constant which should be evaluated once, ahead of time.
=={{header|Java}}==
{{trans|C}}
<syntaxhighlight lang="Java">
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));
}
}
}
</syntaxhighlight>
{{out}}
<pre>
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)
</pre>
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
{{works with|jq}}
'''Works with gojq, the Go implementation of jq'''
'''Preliminaries'''
<syntaxhighlight lang="jq">
# 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) + .;
</syntaxhighlight>
'''The task'''
<syntaxhighlight lang="jq">
# 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
</syntaxhighlight>
'''Invocation:''' jq -ncr -f cordic.jq
{{output}}
<pre>
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
</pre>
=={{header|Julia}}==
Line 41 ⟶ 668:
"""
function cordic(alpha, iteration = 24)
# Fix for the Wikipedia's MATLAB code bug in
newsgn = isodd(Int(floor(alpha / 2π))) ? 1 : -1
alpha < -π/2 && return newsgn * cordic(alpha + π, iteration)
Line 140 ⟶ 767:
+6.0 -0.27941552 (-0.00000002) +0.96017028 (-0.00000001)
</pre>
=={{header|Lua}}==
{{Trans|ALGOL 68}}
<syntaxhighlight lang="lua">
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
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="perl" line>
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
}
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Phix}}==
{{trans|Wren}}
<!--(phixonline)-->
<syntaxhighlight lang="phix">
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
</syntaxhighlight>
{{out}}
As noted above in the commented out line, and mentioned in the Wren entry, testing for <= flips the sign of sin(0).
<pre>
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)
</pre>
=={{header|Python}}==
{{trans|C}}
<syntaxhighlight lang="python">
# 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()
</syntaxhighlight>
{{out}}
<pre>
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)
</pre>
=={{header|Raku}}==
{{trans|XPL0}}
<syntaxhighlight lang="raku" line># 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
} </syntaxhighlight>
{{out}} Same as XPL0 or you may [https://ato.pxeger.com/run?1=XZLfTsIwFMYTL_cUX5aZ8GcrK4OyBSQh6IXhwsRgAgE0YIYSccNuixKCL-KNF_pSPo2nKwvGNm2633fO2Tmn_fiS86fs8_M7S5eO_3Nyk2QL9K-uzy_7KFk9rBLcx5ttGTvDAPC8RWk6DCMb1oDWiNa4jDOUwG1UaPIad8EYo2N5citmNlxbaW678J8OF2vy6KDF_KYX-Fx4XtBqNHw3dDgCFgjhi2a9EXAu6i7BuoIBYSUIJRH0oML9HbmRNlOGLTJqFLDACjaPUGMFBVzmdnWOr4-rdQiqvQuuzHfFr6xBtXpUO6BCJtZg1i7kHpyz_7BEPdItskZwQK1TasUaq86hegSj3GOvNhmmmYzy7iYvMqUQJJOtNSa_srE3jGS-hdmLHigTGofr0qMfJ6soPHxcSBlLs20sY4mOE8AFZ02I7qEmdRsyTLJ1SvlNarMiEpsvkjydjVxF6RLmKeqMLwHzPT_66og3eO_mNDKpkjsbOo4Ndh8neqdyNTT2-nUdHlnx2H4B Attempt This Online!]
=={{header|RPL}}==
Line 195 ⟶ 1,192:
2: "6 constants in memory."
1: { -.91113026188 1 7.07372016661E-2 .960170286655 }
</pre>
=={{header|Rust}}==
{{trans|C}}
<syntaxhighlight lang="Rust">
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());
}
}
</syntaxhighlight>
{{out}}
<pre>
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)
</pre>
=={{header|Sidef}}==
{{trans|Perl}}
<syntaxhighlight lang="ruby">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)
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Wren}}==
{{libheader|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.
<syntaxhighlight lang="wren">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)
}</syntaxhighlight>
{{out}}
<pre>
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)
</pre>
=={{header|XPL0}}==
<syntaxhighlight lang "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);
];
]</syntaxhighlight>
{{out}}
<pre>
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
</pre>
| |||