12
12
13
13
twin_fixed_t twin_sin (twin_angle_t a )
14
14
{
15
- twin_fixed_t sin ;
16
-
17
- /* limit to [0..360) */
18
- a = a & (TWIN_ANGLE_360 - 1 );
19
- /* special case for 90 degrees - no room in table */
20
- if ((a & ~(TWIN_ANGLE_180 )) == TWIN_ANGLE_90 )
21
- sin = TWIN_FIXED_ONE ;
22
- else {
23
- /* mirror second and third quadrant values across y axis */
24
- if (a & TWIN_ANGLE_90 )
25
- a = TWIN_ANGLE_180 - a ;
26
- /* 5th order polynomial approximation */
27
- const uint64_t A1 = 3370945099 , B1 = 2746362156 , C1 = 2339369 ;
28
- const uint64_t A = 16 , n = 10 , p = 32 , q = 31 , r = 3 ;
29
- uint64_t x = a & (TWIN_ANGLE_90 - 1 );
30
- uint64_t y = (C1 * x ) >> n ;
31
- y = B1 - ((x * y ) >> r );
32
- y = x * (y >> n );
33
- y = x * (y >> n );
34
- y = A1 - (y >> (p - q ));
35
- y = x * (y >> n );
36
- y = (y + (1UL << (q - A - 1 ))) >> (q - A ); // Rounding
37
- sin = y ;
38
- }
39
- /* mirror third and fourth quadrant values across x axis */
40
- if (a & TWIN_ANGLE_180 )
41
- sin = - sin ;
42
- return sin ;
15
+ twin_fixed_t sin_val = 0 ;
16
+ twin_sincos (a , & sin_val , NULL );
17
+ return sin_val ;
43
18
}
44
19
45
20
twin_fixed_t twin_cos (twin_angle_t a )
46
21
{
47
- return twin_sin (a + TWIN_ANGLE_90 );
22
+ twin_fixed_t cos_val = 0 ;
23
+ twin_sincos (a , NULL , & cos_val );
24
+ return cos_val ;
48
25
}
49
26
50
27
twin_fixed_t twin_tan (twin_angle_t a )
@@ -63,9 +40,63 @@ twin_fixed_t twin_tan(twin_angle_t a)
63
40
return ((s << 15 ) / c ) << 1 ;
64
41
}
65
42
66
- void twin_sincos ( twin_angle_t a , twin_fixed_t * sin , twin_fixed_t * cos )
43
+ static inline twin_fixed_t sin_poly ( twin_angle_t x )
67
44
{
68
- * sin = twin_sin (a );
69
- * cos = twin_cos (a );
45
+ /* S(x) = x * 2^(-n) * (A1 - 2 ^ (q-p) * x * (2^-n) * x * 2^(-n) * (B1 - 2 ^
46
+ * (-r) * x * 2 ^ (-n) * C1 * x)) * 2 ^ (a-q)
47
+ * @n: the angle scale
48
+ * @A: the amplitude
49
+ * @p,q,r: the scaling factor
50
+ *
51
+ * A1 = 2^q * a5, B1 = 2 ^ p * b5, C1 = 2 ^ (r+p-n) * c5
52
+ * where a5, b5, c5 are the coefficients for 5th-order polynomial
53
+ * a5 = 4 * (3 / pi - 9 / 16)
54
+ * b5 = 2 * a5 - 5 / 2
55
+ * c5 = a5 - 3 / 2
56
+ */
57
+ const uint64_t A = 16 , n = 10 , p = 32 , q = 31 , r = 3 ;
58
+ const uint64_t A1 = 3370945099 , B1 = 2746362156 , C1 = 2339369 ;
59
+ uint64_t y = (C1 * x ) >> n ;
60
+ y = B1 - ((x * y ) >> r );
61
+ y = x * (y >> n );
62
+ y = x * (y >> n );
63
+ y = A1 - (y >> (p - q ));
64
+ y = x * (y >> n );
65
+ y = (y + (1UL << (q - A - 1 ))) >> (q - A ); // Rounding
66
+ return y ;
70
67
}
71
68
69
+ void twin_sincos (twin_angle_t a , twin_fixed_t * sin , twin_fixed_t * cos )
70
+ {
71
+ twin_fixed_t sin_val = 0 , cos_val = 0 ;
72
+
73
+ /* limit to [0..360) */
74
+ a = a & (TWIN_ANGLE_360 - 1 );
75
+ int c = a > TWIN_ANGLE_90 && a < TWIN_ANGLE_270 ;
76
+ /* special case for 90 degrees - no room in table */
77
+ if ((a & ~(TWIN_ANGLE_180 )) == TWIN_ANGLE_90 ) {
78
+ sin_val = TWIN_FIXED_ONE ;
79
+ cos_val = 0 ;
80
+ } else {
81
+ /* mirror second and third quadrant values across y axis */
82
+ if (a & TWIN_ANGLE_90 )
83
+ a = TWIN_ANGLE_180 - a ;
84
+ twin_angle_t x = a & (TWIN_ANGLE_90 - 1 );
85
+ if (sin )
86
+ sin_val = sin_poly (x );
87
+ if (cos )
88
+ cos_val = sin_poly (TWIN_ANGLE_90 - x );
89
+ }
90
+ if (sin ) {
91
+ /* mirror third and fourth quadrant values across x axis */
92
+ if (a & TWIN_ANGLE_180 )
93
+ sin_val = - sin_val ;
94
+ * sin = sin_val ;
95
+ }
96
+ if (cos ) {
97
+ /* mirror first and fourth quadrant values across y axis */
98
+ if (c )
99
+ cos_val = - cos_val ;
100
+ * cos = cos_val ;
101
+ }
102
+ }
0 commit comments