1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #ifdef __RESTRICT
31 #define restrict _Restrict
32 #else
33 #define restrict
34 #endif
35
36 void
37 __vatanf( int n, float * restrict x, int stridex, float * restrict y, int stridey )
38 {
39 extern const double __vlibm_TBL_atan1[];
40 double conup0, conup1, conup2;
41 float dummy, ansf;
42 float f0, f1, f2;
43 float ans0, ans1, ans2;
44 float poly0, poly1, poly2;
45 float sign0, sign1, sign2;
46 int intf, intz, argcount;
47 int index0, index1, index2;
48 float z,*yaddr0,*yaddr1,*yaddr2;
49 int *pz = (int *) &z;
50 #ifdef UNROLL4
51 double conup3;
52 int index3;
53 float f3, ans3, poly3, sign3, *yaddr3;
54 #endif
55
56 /* Power series atan(x) = x + p1*x**3 + p2*x**5 + p3*x**7
57 * Error = -3.08254E-18 On the interval |x| < 1/64 */
58
59 static const float p1 = -0.33329644f /* -3.333333333329292858E-01f */ ;
60 static const float pone = 1.0f;
61
62 if( n <= 0 ) return; /* if no. of elements is 0 or neg, do nothing */
63 do
64 {
65 LOOP0:
66
67 intf = *(int *) x; /* upper half of x, as integer */
68 f0 = *x;
69 sign0 = pone;
70 if (intf < 0) {
71 intf = intf & ~0x80000000; /* abs(upper argument) */
72 f0 = -f0;
73 sign0 = -sign0;
74 }
75
76 if( (intf > 0x5B000000) || (intf < 0x31800000) ) /* filter out special cases */
77 {
78 if( intf > 0x7f800000 )
79 {
80 ansf = f0- f0; /* return NaN if x=NaN*/
81 }
82 else if( intf < 0x31800000 ) /* avoid underflow for small arg */
83 {
84 dummy = 1.0e37 + f0;
85 dummy = dummy;
86 ansf = f0;
87 }
88 else if( intf > 0x5B000000 ) /* avoid underflow for big arg */
89 {
90 index0= 2;
91 ansf = __vlibm_TBL_atan1[index0];/* pi/2 up */
92 }
93 *y = sign0*ansf; /* store answer, with sign bit */
94 x += stridex;
95 y += stridey;
96 argcount = 0; /* initialize argcount */
97 if ( --n <=0 ) break; /* we are done */
98 goto LOOP0; /* otherwise, examine next arg */
99 }
100
101 if (intf > 0x42800000) /* if(|x| > 64 */
102 {
103 f0 = -pone/f0;
104 index0 = 2; /* point to pi/2 upper, lower */
105 }
106 else if( intf >= 0x3C800000 ) /* if |x| >= (1/64)... */
107 {
108 intz = (intf + 0x00040000) & 0x7ff80000;/* round arg, keep upper */
109 pz[0] = intz; /* store as a float (z) */
110 f0 = (f0 - z)/(pone + f0*z);
111 index0 = (intz - 0x3C800000) >> 18; /* (index >> 19) << 1) */
112 index0 = index0+ 4; /* skip over 0,0,pi/2,pi/2 */
113 }
114 else /* |x| < 1/64 */
115 {
116 index0 = 0; /* points to 0,0 in table */
117 }
118 yaddr0 = y; /* address to store this answer */
119 x += stridex; /* point to next arg */
120 y += stridey; /* point to next result */
121 argcount = 1; /* we now have 1 good argument */
122 if ( --n <=0 )
123 {
124 goto UNROLL; /* finish up with 1 good arg */
125 }
126
127 /*--------------------------------------------------------------------------*/
128 /*--------------------------------------------------------------------------*/
129 /*--------------------------------------------------------------------------*/
130
131 LOOP1:
132
133 intf = *(int *) x; /* upper half of x, as integer */
134 f1 = *x;
135 sign1 = pone;
136 if (intf < 0) {
137 intf = intf & ~0x80000000; /* abs(upper argument) */
138 f1 = -f1;
139 sign1 = -sign1;
140 }
141
142 if( (intf > 0x5B000000) || (intf < 0x31800000) ) /* filter out special cases */
143 {
144 if( intf > 0x7f800000 )
145 {
146 ansf = f1 - f1; /* return NaN if x=NaN*/
147 }
148 else if( intf < 0x31800000 ) /* avoid underflow for small arg */
149 {
150 dummy = 1.0e37 + f1;
151 dummy = dummy;
152 ansf = f1;
153 }
154 else if( intf > 0x5B000000 ) /* avoid underflow for big arg */
155 {
156 index1 = 2;
157 ansf = __vlibm_TBL_atan1[index1] ;/* pi/2 up */
158 }
159 *y = sign1 * ansf; /* store answer, with sign bit */
160 x += stridex;
161 y += stridey;
162 argcount = 1; /* we still have 1 good arg */
163 if ( --n <=0 )
164 {
165 goto UNROLL; /* finish up with 1 good arg */
166 }
167 goto LOOP1; /* otherwise, examine next arg */
168 }
169
170 if (intf > 0x42800000) /* if(|x| > 64 */
171 {
172 f1 = -pone/f1;
173 index1 = 2; /* point to pi/2 upper, lower */
174 }
175 else if( intf >= 0x3C800000 ) /* if |x| >= (1/64)... */
176 {
177 intz = (intf + 0x00040000) & 0x7ff80000;/* round arg, keep upper */
178 pz[0] = intz; /* store as a float (z) */
179 f1 = (f1 - z)/(pone + f1*z);
180 index1 = (intz - 0x3C800000) >> 18; /* (index >> 19) << 1) */
181 index1 = index1 + 4; /* skip over 0,0,pi/2,pi/2 */
182 }
183 else
184 {
185 index1 = 0; /* points to 0,0 in table */
186 }
187
188 yaddr1 = y; /* address to store this answer */
189 x += stridex; /* point to next arg */
190 y += stridey; /* point to next result */
191 argcount = 2; /* we now have 2 good arguments */
192 if ( --n <=0 )
193 {
194 goto UNROLL; /* finish up with 2 good args */
195 }
196
197 /*--------------------------------------------------------------------------*/
198 /*--------------------------------------------------------------------------*/
199 /*--------------------------------------------------------------------------*/
200
201 LOOP2:
202
203 intf = *(int *) x; /* upper half of x, as integer */
204 f2 = *x;
205 sign2 = pone;
206 if (intf < 0) {
207 intf = intf & ~0x80000000; /* abs(upper argument) */
208 f2 = -f2;
209 sign2 = -sign2;
210 }
211
212 if( (intf > 0x5B000000) || (intf < 0x31800000) ) /* filter out special cases */
213 {
214 if( intf > 0x7f800000 )
215 {
216 ansf = f2 - f2; /* return NaN if x=NaN*/
217 }
218 else if( intf < 0x31800000 ) /* avoid underflow for small arg */
219 {
220 dummy = 1.0e37 + f2;
221 dummy = dummy;
222 ansf = f2;
223 }
224 else if( intf > 0x5B000000 ) /* avoid underflow for big arg */
225 {
226 index2 = 2;
227 ansf = __vlibm_TBL_atan1[index2] ;/* pi/2 up */
228 }
229 *y = sign2 * ansf; /* store answer, with sign bit */
230 x += stridex;
231 y += stridey;
232 argcount = 2; /* we still have 2 good args */
233 if ( --n <=0 )
234 {
235 goto UNROLL; /* finish up with 2 good args */
236 }
237 goto LOOP2; /* otherwise, examine next arg */
238 }
239
240 if (intf > 0x42800000) /* if(|x| > 64 */
241 {
242 f2 = -pone/f2;
243 index2 = 2; /* point to pi/2 upper, lower */
244 }
245 else if( intf >= 0x3C800000 ) /* if |x| >= (1/64)... */
246 {
247 intz = (intf + 0x00040000) & 0x7ff80000;/* round arg, keep upper */
248 pz[0] = intz; /* store as a float (z) */
249 f2 = (f2 - z)/(pone + f2*z);
250 index2 = (intz - 0x3C800000) >> 18; /* (index >> 19) << 1) */
251 index2 = index2 + 4; /* skip over 0,0,pi/2,pi/2 */
252 }
253 else
254 {
255 index2 = 0; /* points to 0,0 in table */
256 }
257 yaddr2 = y; /* address to store this answer */
258 x += stridex; /* point to next arg */
259 y += stridey; /* point to next result */
260 argcount = 3; /* we now have 3 good arguments */
261 if ( --n <=0 )
262 {
263 goto UNROLL; /* finish up with 2 good args */
264 }
265
266
267 /*--------------------------------------------------------------------------*/
268 /*--------------------------------------------------------------------------*/
269 /*--------------------------------------------------------------------------*/
270
271 #ifdef UNROLL4
272 LOOP3:
273
274 intf = *(int *) x; /* upper half of x, as integer */
275 f3 = *x;
276 sign3 = pone;
277 if (intf < 0) {
278 intf = intf & ~0x80000000; /* abs(upper argument) */
279 f3 = -f3;
280 sign3 = -sign3;
281 }
282
283 if( (intf > 0x5B000000) || (intf < 0x31800000) ) /* filter out special cases */
284 {
285 if( intf > 0x7f800000 )
286 {
287 ansf = f3 - f3; /* return NaN if x=NaN*/
288 }
289 else if( intf < 0x31800000 ) /* avoid underflow for small arg */
290 {
291 dummy = 1.0e37 + f3;
292 dummy = dummy;
293 ansf = f3;
294 }
295 else if( intf > 0x5B000000 ) /* avoid underflow for big arg */
296 {
297 index3 = 2;
298 ansf = __vlibm_TBL_atan1[index3] ;/* pi/2 up */
299 }
300 *y = sign3 * ansf; /* store answer, with sign bit */
301 x += stridex;
302 y += stridey;
303 argcount = 3; /* we still have 3 good args */
304 if ( --n <=0 )
305 {
306 goto UNROLL; /* finish up with 3 good args */
307 }
308 goto LOOP3; /* otherwise, examine next arg */
309 }
310
311 if (intf > 0x42800000) /* if(|x| > 64 */
312 {
313 n3 = -pone;
314 d3 = f3;
315 f3 = n3/d3;
316 index3 = 2; /* point to pi/2 upper, lower */
317 }
318 else if( intf >= 0x3C800000 ) /* if |x| >= (1/64)... */
319 {
320 intz = (intf + 0x00040000) & 0x7ff80000;/* round arg, keep upper */
321 pz[0] = intz; /* store as a float (z) */
322 n3 = (f3 - z);
323 d3 = (pone + f3*z); /* get reduced argument */
324 f3 = n3/d3;
325 index3 = (intz - 0x3C800000) >> 18; /* (index >> 19) << 1) */
326 index3 = index3 + 4; /* skip over 0,0,pi/2,pi/2 */
327 }
328 else
329 {
330 n3 = f3;
331 d3 = pone;
332 index3 = 0; /* points to 0,0 in table */
333 }
334 yaddr3 = y; /* address to store this answer */
335 x += stridex; /* point to next arg */
336 y += stridey; /* point to next result */
337 argcount = 4; /* we now have 4 good arguments */
338 if ( --n <=0 )
339 {
340 goto UNROLL; /* finish up with 3 good args */
341 }
342 #endif /* UNROLL4 */
343
344 /* here is the n-way unrolled section,
345 but we may actually have less than n
346 arguments at this point
347 */
348
349 UNROLL:
350
351 #ifdef UNROLL4
352 if (argcount == 4)
353 {
354 conup0 = __vlibm_TBL_atan1[index0];
355 conup1 = __vlibm_TBL_atan1[index1];
356 conup2 = __vlibm_TBL_atan1[index2];
357 conup3 = __vlibm_TBL_atan1[index3];
358 poly0 = p1*f0*f0*f0 + f0;
359 ans0 = sign0 * (float)(conup0 + poly0);
360 poly1 = p1*f1*f1*f1 + f1;
361 ans1 = sign1 * (float)(conup1 + poly1);
362 poly2 = p1*f2*f2*f2 + f2;
363 ans2 = sign2 * (float)(conup2 + poly2);
364 poly3 = p1*f3*f3*f3 + f3;
365 ans3 = sign3 * (float)(conup3 + poly3);
366 *yaddr0 = ans0;
367 *yaddr1 = ans1;
368 *yaddr2 = ans2;
369 *yaddr3 = ans3;
370 }
371 else
372 #endif
373 if (argcount == 3)
374 {
375 conup0 = __vlibm_TBL_atan1[index0];
376 conup1 = __vlibm_TBL_atan1[index1];
377 conup2 = __vlibm_TBL_atan1[index2];
378 poly0 = p1*f0*f0*f0 + f0;
379 poly1 = p1*f1*f1*f1 + f1;
380 poly2 = p1*f2*f2*f2 + f2;
381 ans0 = sign0 * (float)(conup0 + poly0);
382 ans1 = sign1 * (float)(conup1 + poly1);
383 ans2 = sign2 * (float)(conup2 + poly2);
384 *yaddr0 = ans0;
385 *yaddr1 = ans1;
386 *yaddr2 = ans2;
387 }
388 else
389 if (argcount == 2)
390 {
391 conup0 = __vlibm_TBL_atan1[index0];
392 conup1 = __vlibm_TBL_atan1[index1];
393 poly0 = p1*f0*f0*f0 + f0;
394 poly1 = p1*f1*f1*f1 + f1;
395 ans0 = sign0 * (float)(conup0 + poly0);
396 ans1 = sign1 * (float)(conup1 + poly1);
397 *yaddr0 = ans0;
398 *yaddr1 = ans1;
399 }
400 else
401 if (argcount == 1)
402 {
403 conup0 = __vlibm_TBL_atan1[index0];
404 poly0 = p1*f0*f0*f0 + f0;
405 ans0 = sign0 * (float)(conup0 + poly0);
406 *yaddr0 = ans0;
407 }
408
409 } while (n > 0);
410
411 }