?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
2.9
[glibc/nacl-glibc.git] / sysdeps / ieee754 / dbl-64 / s_tan.c
blobcf8d4d02670f8143b153bf54e5cc825dc2da53e5
1 /*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 Free Software Foundation
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19 */
20 /*********************************************************************/
21 /* MODULE_NAME: utan.c */
22 /* */
23 /* FUNCTIONS: utan */
24 /* tanMp */
25 /* */
26 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h utan.h */
27 /* branred.c sincos32.c mptan.c */
28 /* utan.tbl */
29 /* */
30 /* An ultimate tan routine. Given an IEEE double machine number x */
31 /* it computes the correctly rounded (to nearest) value of tan(x). */
32 /* Assumption: Machine arithmetic operations are performed in */
33 /* round to nearest mode of IEEE 754 standard. */
34 /* */
35 /*********************************************************************/
36 #include "endian.h"
37 #include "dla.h"
38 #include "mpa.h"
39 #include "MathLib.h"
40 #include "math.h"
41
42 static double tanMp(double);
43 void __mptan(double, mp_no *, int);
44
45 double tan(double x) {
46 #include "utan.h"
47 #include "utan.tbl"
48
49 int ux,i,n;
50 double a,da,a2,b,db,c,dc,c1,cc1,c2,cc2,c3,cc3,fi,ffi,gi,pz,s,sy,
51 t,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2;
52 int p;
53 number num,v;
54 mp_no mpa,mpt1,mpt2;
55 #if 0
56 mp_no mpy;
57 #endif
58
59 int __branred(double, double *, double *);
60 int __mpranred(double, mp_no *, int);
61
62 /* x=+-INF, x=NaN */
63 num.d = x; ux = num.i[HIGH_HALF];
64 if ((ux&0x7ff00000)==0x7ff00000) return x-x;
65
66 w=(x<ZERO) ? -x : x;
67
68 /* (I) The case abs(x) <= 1.259e-8 */
69 if (w<=g1.d) return x;
70
71 /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
72 if (w<=g2.d) {
73
74 /* First stage */
75 x2 = x*x;
76 t2 = x*x2*(d3.d+x2*(d5.d+x2*(d7.d+x2*(d9.d+x2*d11.d))));
77 if ((y=x+(t2-u1.d*t2)) == x+(t2+u1.d*t2)) return y;
78
79 /* Second stage */
80 c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
81 x2*a27.d))))));
82 EMULV(x,x,x2,xx2,t1,t2,t3,t4,t5)
83 ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
84 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
85 ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
86 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
87 ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
88 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
89 ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
90 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
91 ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
92 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
93 ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
94 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
95 MUL2(x ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
96 ADD2(x ,zero.d,c2,cc2,c1,cc1,t1,t2)
97 if ((y=c1+(cc1-u2.d*c1)) == c1+(cc1+u2.d*c1)) return y;
98 return tanMp(x);
99 }
100
101 /* (III) The case 0.0608 < abs(x) <= 0.787 */
102 if (w<=g3.d) {
103
104 /* First stage */
105 i = ((int) (mfftnhf.d+TWO8*w));
106 z = w-xfg[i][0].d; z2 = z*z; s = (x<ZERO) ? MONE : ONE;
107 pz = z+z*z2*(e0.d+z2*e1.d);
108 fi = xfg[i][1].d; gi = xfg[i][2].d; t2 = pz*(gi+fi)/(gi-pz);
109 if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d)) return (s*y);
110 t3 = (t2<ZERO) ? -t2 : t2;
111 if ((y=fi+(t2-(t4=fi*ua3.d+t3*ub3.d)))==fi+(t2+t4)) return (s*y);
112
113 /* Second stage */
114 ffi = xfg[i][3].d;
115 c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
116 EMULV(z,z,z2,zz2,t1,t2,t3,t4,t5)
117 ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
118 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
119 ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
120 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
121 MUL2(z ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
122 ADD2(z ,zero.d,c2,cc2,c1,cc1,t1,t2)
123
124 ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
125 MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
126 SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
127 DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
128
129 if ((y=c3+(cc3-u4.d*c3))==c3+(cc3+u4.d*c3)) return (s*y);
130 return tanMp(x);
131 }
132
133 /* (---) The case 0.787 < abs(x) <= 25 */
134 if (w<=g4.d) {
135 /* Range reduction by algorithm i */
136 t = (x*hpinv.d + toint.d);
137 xn = t - toint.d;
138 v.d = t;
139 t1 = (x - xn*mp1.d) - xn*mp2.d;
140 n =v.i[LOW_HALF] & 0x00000001;
141 da = xn*mp3.d;
142 a=t1-da;
143 da = (t1-a)-da;
144 if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
145 else {ya= a; yya= da; sy= ONE;}
146
147 /* (IV),(V) The case 0.787 < abs(x) <= 25, abs(y) <= 1e-7 */
148 if (ya<=gy1.d) return tanMp(x);
149
150 /* (VI) The case 0.787 < abs(x) <= 25, 1e-7 < abs(y) <= 0.0608 */
151 if (ya<=gy2.d) {
152 a2 = a*a;
153 t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
154 if (n) {
155 /* First stage -cot */
156 EADD(a,t2,b,db)
157 DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
158 if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c)) return (-y); }
159 else {
160 /* First stage tan */
161 if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a)) return y; }
162 /* Second stage */
163 /* Range reduction by algorithm ii */
164 t = (x*hpinv.d + toint.d);
165 xn = t - toint.d;
166 v.d = t;
167 t1 = (x - xn*mp1.d) - xn*mp2.d;
168 n =v.i[LOW_HALF] & 0x00000001;
169 da = xn*pp3.d;
170 t=t1-da;
171 da = (t1-t)-da;
172 t1 = xn*pp4.d;
173 a = t - t1;
174 da = ((t-a)-t1)+da;
175
176 /* Second stage */
177 EADD(a,da,t1,t2) a=t1; da=t2;
178 MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
179 c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
180 x2*a27.d))))));
181 ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
182 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
183 ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
184 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
185 ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
186 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
187 ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
188 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
189 ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
190 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
191 ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
192 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
193 MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
194 ADD2(a ,da ,c2,cc2,c1,cc1,t1,t2)
195
196 if (n) {
197 /* Second stage -cot */
198 DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
199 if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2)) return (-y); }
200 else {
201 /* Second stage tan */
202 if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1)) return y; }
203 return tanMp(x);
204 }
205
206 /* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */
207
208 /* First stage */
209 i = ((int) (mfftnhf.d+TWO8*ya));
210 z = (z0=(ya-xfg[i][0].d))+yya; z2 = z*z;
211 pz = z+z*z2*(e0.d+z2*e1.d);
212 fi = xfg[i][1].d; gi = xfg[i][2].d;
213
214 if (n) {
215 /* -cot */
216 t2 = pz*(fi+gi)/(fi+pz);
217 if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d)) return (-sy*y);
218 t3 = (t2<ZERO) ? -t2 : t2;
219 if ((y=gi-(t2-(t4=gi*ua10.d+t3*ub10.d)))==gi-(t2+t4)) return (-sy*y); }
220 else {
221 /* tan */
222 t2 = pz*(gi+fi)/(gi-pz);
223 if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d)) return (sy*y);
224 t3 = (t2<ZERO) ? -t2 : t2;
225 if ((y=fi+(t2-(t4=fi*ua9.d+t3*ub9.d)))==fi+(t2+t4)) return (sy*y); }
226
227 /* Second stage */
228 ffi = xfg[i][3].d;
229 EADD(z0,yya,z,zz)
230 MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
231 c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
232 ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
233 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
234 ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
235 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
236 MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
237 ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
238
239 ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
240 MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
241 SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
242
243 if (n) {
244 /* -cot */
245 DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
246 if ((y=c3+(cc3-u12.d*c3))==c3+(cc3+u12.d*c3)) return (-sy*y); }
247 else {
248 /* tan */
249 DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
250 if ((y=c3+(cc3-u11.d*c3))==c3+(cc3+u11.d*c3)) return (sy*y); }
251
252 return tanMp(x);
253 }
254
255 /* (---) The case 25 < abs(x) <= 1e8 */
256 if (w<=g5.d) {
257 /* Range reduction by algorithm ii */
258 t = (x*hpinv.d + toint.d);
259 xn = t - toint.d;
260 v.d = t;
261 t1 = (x - xn*mp1.d) - xn*mp2.d;
262 n =v.i[LOW_HALF] & 0x00000001;
263 da = xn*pp3.d;
264 t=t1-da;
265 da = (t1-t)-da;
266 t1 = xn*pp4.d;
267 a = t - t1;
268 da = ((t-a)-t1)+da;
269 EADD(a,da,t1,t2) a=t1; da=t2;
270 if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
271 else {ya= a; yya= da; sy= ONE;}
272
273 /* (+++) The case 25 < abs(x) <= 1e8, abs(y) <= 1e-7 */
274 if (ya<=gy1.d) return tanMp(x);
275
276 /* (VIII) The case 25 < abs(x) <= 1e8, 1e-7 < abs(y) <= 0.0608 */
277 if (ya<=gy2.d) {
278 a2 = a*a;
279 t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
280 if (n) {
281 /* First stage -cot */
282 EADD(a,t2,b,db)
283 DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
284 if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c)) return (-y); }
285 else {
286 /* First stage tan */
287 if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a)) return y; }
288
289 /* Second stage */
290 MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
291 c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
292 x2*a27.d))))));
293 ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
294 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
295 ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
296 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
297 ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
298 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
299 ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
300 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
301 ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
302 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
303 ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
304 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
305 MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
306 ADD2(a ,da ,c2,cc2,c1,cc1,t1,t2)
307
308 if (n) {
309 /* Second stage -cot */
310 DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
311 if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2)) return (-y); }
312 else {
313 /* Second stage tan */
314 if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1)) return (y); }
315 return tanMp(x);
316 }
317
318 /* (IX) The case 25 < abs(x) <= 1e8, 0.0608 < abs(y) <= 0.787 */
319 /* First stage */
320 i = ((int) (mfftnhf.d+TWO8*ya));
321 z = (z0=(ya-xfg[i][0].d))+yya; z2 = z*z;
322 pz = z+z*z2*(e0.d+z2*e1.d);
323 fi = xfg[i][1].d; gi = xfg[i][2].d;
324
325 if (n) {
326 /* -cot */
327 t2 = pz*(fi+gi)/(fi+pz);
328 if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d)) return (-sy*y);
329 t3 = (t2<ZERO) ? -t2 : t2;
330 if ((y=gi-(t2-(t4=gi*ua18.d+t3*ub18.d)))==gi-(t2+t4)) return (-sy*y); }
331 else {
332 /* tan */
333 t2 = pz*(gi+fi)/(gi-pz);
334 if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d)) return (sy*y);
335 t3 = (t2<ZERO) ? -t2 : t2;
336 if ((y=fi+(t2-(t4=fi*ua17.d+t3*ub17.d)))==fi+(t2+t4)) return (sy*y); }
337
338 /* Second stage */
339 ffi = xfg[i][3].d;
340 EADD(z0,yya,z,zz)
341 MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
342 c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
343 ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
344 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
345 ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
346 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
347 MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
348 ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
349
350 ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
351 MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
352 SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
353
354 if (n) {
355 /* -cot */
356 DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
357 if ((y=c3+(cc3-u20.d*c3))==c3+(cc3+u20.d*c3)) return (-sy*y); }
358 else {
359 /* tan */
360 DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
361 if ((y=c3+(cc3-u19.d*c3))==c3+(cc3+u19.d*c3)) return (sy*y); }
362 return tanMp(x);
363 }
364
365 /* (---) The case 1e8 < abs(x) < 2**1024 */
366 /* Range reduction by algorithm iii */
367 n = (__branred(x,&a,&da)) & 0x00000001;
368 EADD(a,da,t1,t2) a=t1; da=t2;
369 if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
370 else {ya= a; yya= da; sy= ONE;}
371
372 /* (+++) The case 1e8 < abs(x) < 2**1024, abs(y) <= 1e-7 */
373 if (ya<=gy1.d) return tanMp(x);
374
375 /* (X) The case 1e8 < abs(x) < 2**1024, 1e-7 < abs(y) <= 0.0608 */
376 if (ya<=gy2.d) {
377 a2 = a*a;
378 t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
379 if (n) {
380 /* First stage -cot */
381 EADD(a,t2,b,db)
382 DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
383 if ((y=c+(dc-u22.d*c))==c+(dc+u22.d*c)) return (-y); }
384 else {
385 /* First stage tan */
386 if ((y=a+(t2-u21.d*a))==a+(t2+u21.d*a)) return y; }
387
388 /* Second stage */
389 /* Reduction by algorithm iv */
390 p=10; n = (__mpranred(x,&mpa,p)) & 0x00000001;
391 __mp_dbl(&mpa,&a,p); __dbl_mp(a,&mpt1,p);
392 __sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p);
393
394 MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
395 c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
396 x2*a27.d))))));
397 ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
398 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
399 ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
400 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
401 ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
402 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
403 ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
404 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
405 ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
406 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
407 ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
408 MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
409 MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
410 ADD2(a ,da ,c2,cc2,c1,cc1,t1,t2)
411
412 if (n) {
413 /* Second stage -cot */
414 DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
415 if ((y=c2+(cc2-u24.d*c2)) == c2+(cc2+u24.d*c2)) return (-y); }
416 else {
417 /* Second stage tan */
418 if ((y=c1+(cc1-u23.d*c1)) == c1+(cc1+u23.d*c1)) return y; }
419 return tanMp(x);
420 }
421
422 /* (XI) The case 1e8 < abs(x) < 2**1024, 0.0608 < abs(y) <= 0.787 */
423 /* First stage */
424 i = ((int) (mfftnhf.d+TWO8*ya));
425 z = (z0=(ya-xfg[i][0].d))+yya; z2 = z*z;
426 pz = z+z*z2*(e0.d+z2*e1.d);
427 fi = xfg[i][1].d; gi = xfg[i][2].d;
428
429 if (n) {
430 /* -cot */
431 t2 = pz*(fi+gi)/(fi+pz);
432 if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d)) return (-sy*y);
433 t3 = (t2<ZERO) ? -t2 : t2;
434 if ((y=gi-(t2-(t4=gi*ua26.d+t3*ub26.d)))==gi-(t2+t4)) return (-sy*y); }
435 else {
436 /* tan */
437 t2 = pz*(gi+fi)/(gi-pz);
438 if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d)) return (sy*y);
439 t3 = (t2<ZERO) ? -t2 : t2;
440 if ((y=fi+(t2-(t4=fi*ua25.d+t3*ub25.d)))==fi+(t2+t4)) return (sy*y); }
441
442 /* Second stage */
443 ffi = xfg[i][3].d;
444 EADD(z0,yya,z,zz)
445 MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
446 c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
447 ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
448 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
449 ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
450 MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
451 MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
452 ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
453
454 ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
455 MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
456 SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
457
458 if (n) {
459 /* -cot */
460 DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
461 if ((y=c3+(cc3-u28.d*c3))==c3+(cc3+u28.d*c3)) return (-sy*y); }
462 else {
463 /* tan */
464 DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
465 if ((y=c3+(cc3-u27.d*c3))==c3+(cc3+u27.d*c3)) return (sy*y); }
466 return tanMp(x);
467 }
468
469
470 /* multiple precision stage */
471 /* Convert x to multi precision number,compute tan(x) by mptan() routine */
472 /* and converts result back to double */
473 static double tanMp(double x)
474 {
475 int p;
476 double y;
477 mp_no mpy;
478 p=32;
479 __mptan(x, &mpy, p);
480 __mp_dbl(&mpy,&y,p);
481 return y;
482 }
483
484 #ifdef NO_LONG_DOUBLE
485 weak_alias (tan, tanl)
486 #endif
Port of glibc to Google Native Client (NaCl)