-
Notifications
You must be signed in to change notification settings - Fork 21
/
mysincosf.h
executable file
·82 lines (70 loc) · 1.55 KB
/
mysincosf.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
/*
* Project : SIMD_Utils
* Version : 0.2.5
* Author : JishinMaster
* Licence : BSD-2
*/
#pragma once
static inline int mysincosf(float xx, float *s, float *c)
{
float x, y, y1, y2, z;
int j, sign_sin, sign_cos;
sign_sin = 1;
sign_cos = 1;
x = xx;
if (xx < 0) {
sign_sin = -1;
x = -xx;
}
if (x > T24M1) {
return (0.0f);
}
z = FOPI * x; /* integer part of x/(PI/4) */
j = (int)z;
y = (float)j;
/* map zeros to origin */
if (j & 1) {
j += 1;
y += 1.0f;
}
j &= 7; /* octant modulo 360 degrees */
/* reflect in x axis */
if (j > 3) {
sign_sin = -sign_sin;
sign_cos = -sign_cos;
j -= 4;
}
if (j > 1)
sign_cos = -sign_cos;
if (x > lossth) {
x = x - y * PIO4F;
} else {
/* Extended precision modular arithmetic */
x = ((x + y * minus_cephes_DP1) + y * minus_cephes_DP2) + y * minus_cephes_DP3;
}
/*einits();*/
z = x * x;
/* measured relative error in +/- pi/4 is 7.8e-8 */
y1 = ((coscof[0] * z + coscof[1]) * z + coscof[2]) * z * z;
y1 -= 0.5f * z;
y1 += 1.0f;
/* Theoretical relative error = 3.8e-9 in [-pi/4, +pi/4] */
y2 = ((sincof[0] * z + sincof[1]) * z + sincof[2]) * z * x;
y2 += x;
if ((j == 1) || (j == 2)) {
*s = y1;
*c = y2;
} else {
*s = y2;
*c = y1;
}
// COS
/*einitd();*/
if (sign_sin < 0) {
*s = -(*s);
}
if (sign_cos < 0) {
*c = -(*c);
}
return 0;
}