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: #include "math.h"
48: #include "math_private.h"
49:
50: static const double xxx[] = {
51: 3.33333333333334091986e-01,
52: 1.33333333333201242699e-01,
53: 5.39682539762260521377e-02,
54: 2.18694882948595424599e-02,
55: 8.86323982359930005737e-03,
56: 3.59207910759131235356e-03,
57: 1.45620945432529025516e-03,
58: 5.88041240820264096874e-04,
59: 2.46463134818469906812e-04,
60: 7.81794442939557092300e-05,
61: 7.14072491382608190305e-05,
62: -1.85586374855275456654e-05,
63: 2.59073051863633712884e-05,
64: 1.00000000000000000000e+00,
65: 7.85398163397448278999e-01,
66: 3.06161699786838301793e-17
67: };
68: #define one xxx[13]
69: #define pio4 xxx[14]
70: #define pio4lo xxx[15]
71: #define T xxx
72:
73: double
74: __kernel_tan(double x, double y, int iy)
75: {
76: double z, r, v, w, s;
77: int32_t ix, hx;
78:
79: GET_HIGH_WORD(hx, x);
80: ix = hx & 0x7fffffff;
81: if (ix < 0x3e300000) {
82: if ((int) x == 0) {
83: u_int32_t low;
84: GET_LOW_WORD(low, x);
85: if(((ix | low) | (iy + 1)) == 0)
86: return one / fabs(x);
87: else {
88: if (iy == 1)
89: return x;
90: else {
91: double a, t;
92:
93: z = w = x + y;
94: SET_LOW_WORD(z, 0);
95: v = y - (z - x);
96: t = a = -one / w;
97: SET_LOW_WORD(t, 0);
98: s = one + t * z;
99: return t + a * (s + t * v);
100: }
101: }
102: }
103: }
104: if (ix >= 0x3FE59428) {
105: if (hx < 0) {
106: x = -x;
107: y = -y;
108: }
109: z = pio4 - x;
110: w = pio4lo - y;
111: x = z + w;
112: y = 0.0;
113: }
114: z = x * x;
115: w = z * z;
116: 117: 118: 119: 120:
121: r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] +
122: w * T[11]))));
123: v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] +
124: w * T[12])))));
125: s = z * x;
126: r = y + z * (s * (r + v) + y);
127: r += T[0] * s;
128: w = x + r;
129: if (ix >= 0x3FE59428) {
130: v = (double) iy;
131: return (double) (1 - ((hx >> 30) & 2)) *
132: (v - 2.0 * (x - (w * w / (w + v) - r)));
133: }
134: if (iy == 1)
135: return w;
136: else {
137: 138: 139: 140:
141:
142: double a, t;
143: z = w;
144: SET_LOW_WORD(z, 0);
145: v = r - (z - x);
146: t = a = -1.0 / w;
147: SET_LOW_WORD(t, 0);
148: s = 1.0 + t * z;
149: return t + a * (s + t * v);
150: }
151: }