25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define restrict _Restrict
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define restrict
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__vatan(int n, double * restrict x, int stridex, double * restrict y, int stridey)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis double f, z, ans = 0.0L, ansu, ansl, tmp, poly, conup, conlo, dummy;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis double f1, ans1, ansu1, ansl1, tmp1, poly1, conup1, conlo1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis double f2, ans2, ansu2, ansl2, tmp2, poly2, conup2, conlo2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int index, sign, intf, intflo, intz, argcount;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis extern const double __vlibm_TBL_atan1[];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis extern double fabs(double);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* Power series atan(x) = x + p1*x**3 + p2*x**5 + p3*x**7
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Error = -3.08254E-18 On the interval |x| < 1/64 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* define dummy names for readability. Use parray to help compiler optimize loads */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (n <= 0) return; /* if no. of elements is 0 or neg, do nothing */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intf = HI(x); /* upper half of x, as integer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intflo = LO(x); /* lower half of x, as integer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sign = intf & 0x80000000; /* sign of argument */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intf = intf & ~0x80000000; /* abs(upper argument) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((intf > 0x43600000) || (intf < 0x3e300000)) /* filter out special cases */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ( (intf > 0x7ff00000) || ((intf == 0x7ff00000) && (intflo !=0)))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf < 0x3e300000) /* avoid underflow for small arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf > 0x43600000) /* avoid underflow for big arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ans = __vlibm_TBL_atan1[index] + __vlibm_TBL_atan1[index+1];/* pi/2 up + pi/2 low */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *y = (sign) ? -ans: ans; /* store answer, with sign bit */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (--n <=0) break; /* we are done */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis goto LOOP0; /* otherwise, examine next arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index = 2; /* point to pi/2 upper, lower */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf >= 0x3f900000) /* if |x| >= (1/64)... */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intz = (intf + 0x00008000) & 0x7fff0000;/* round arg, keep upper */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis HI(&z) = intz; /* store as a double (z) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis f = (f - z)/(1.0 + f*z); /* get reduced argument */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index = (intz - 0x3f900000) >> 15; /* (index >> 16) << 1) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index = index + 4; /* skip over 0,0,pi/2,pi/2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis yaddr = y; /* address to store this answer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis argcount = 1; /* we now have 1 good argument */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis f1 = 0.0; /* put dummy values in args 1,2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis goto UNROLL3; /* finish up with 1 good arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /*--------------------------------------------------------------------------*/
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /*--------------------------------------------------------------------------*/
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /*--------------------------------------------------------------------------*/
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intf = HI(x); /* upper half of x, as integer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intflo = LO(x); /* lower half of x, as integer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sign1 = intf & 0x80000000; /* sign of argument */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intf = intf & ~0x80000000; /* abs(upper argument) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((intf > 0x43600000) || (intf < 0x3e300000)) /* filter out special cases */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ( (intf > 0x7ff00000) || ((intf == 0x7ff00000) && (intflo !=0)))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf < 0x3e300000) /* avoid underflow for small arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf > 0x43600000) /* avoid underflow for big arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ans = __vlibm_TBL_atan1[index1] + __vlibm_TBL_atan1[index1+1];/* pi/2 up + pi/2 low */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *y = (sign1) ? -ans: ans; /* store answer, with sign bit */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis argcount = 1; /* we still have 1 good arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis f1 = 0.0; /* put dummy values in args 1,2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis goto UNROLL3; /* finish up with 1 good arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis goto LOOP1; /* otherwise, examine next arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index1 = 2; /* point to pi/2 upper, lower */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf >= 0x3f900000) /* if |x| >= (1/64)... */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intz = (intf + 0x00008000) & 0x7fff0000;/* round arg, keep upper */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis HI(&z) = intz; /* store as a double (z) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis f1 = (f1 - z)/(1.0 + f1*z); /* get reduced argument */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index1 = (intz - 0x3f900000) >> 15; /* (index >> 16) << 1) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index1 = index1 + 4; /* skip over 0,0,pi/2,pi/2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis yaddr1 = y; /* address to store this answer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis argcount = 2; /* we now have 2 good arguments */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis goto UNROLL3; /* finish up with 2 good args */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /*--------------------------------------------------------------------------*/
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /*--------------------------------------------------------------------------*/
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /*--------------------------------------------------------------------------*/
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intf = HI(x); /* upper half of x, as integer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intflo = LO(x); /* lower half of x, as integer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sign2 = intf & 0x80000000; /* sign of argument */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intf = intf & ~0x80000000; /* abs(upper argument) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((intf > 0x43600000) || (intf < 0x3e300000)) /* filter out special cases */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ( (intf > 0x7ff00000) || ((intf == 0x7ff00000) && (intflo !=0)))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf < 0x3e300000) /* avoid underflow for small arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf > 0x43600000) /* avoid underflow for big arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ans = __vlibm_TBL_atan1[index2] + __vlibm_TBL_atan1[index2+1];/* pi/2 up + pi/2 low */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *y = (sign2) ? -ans: ans; /* store answer, with sign bit */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis argcount = 2; /* we still have 2 good args */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis goto UNROLL3; /* finish up with 2 good args */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis goto LOOP2; /* otherwise, examine next arg */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index2 = 2; /* point to pi/2 upper, lower */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (intf >= 0x3f900000) /* if |x| >= (1/64)... */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis intz = (intf + 0x00008000) & 0x7fff0000;/* round arg, keep upper */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis HI(&z) = intz; /* store as a double (z) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis f2 = (f2 - z)/(1.0 + f2*z); /* get reduced argument */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index2 = (intz - 0x3f900000) >> 15; /* (index >> 16) << 1) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis index2 = index2 + 4; /* skip over 0,0,pi/2,pi/2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis yaddr2 = y; /* address to store this answer */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis argcount = 3; /* we now have 3 good arguments */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* here is the 3 way unrolled section,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis note, we may actually only have
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 1,2, or 3 'real' arguments at this point
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis conup = __vlibm_TBL_atan1[index ]; /* upper table */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis conup1 = __vlibm_TBL_atan1[index1]; /* upper table */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis conup2 = __vlibm_TBL_atan1[index2]; /* upper table */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis conlo = __vlibm_TBL_atan1[index +1]; /* lower table */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis conlo1 = __vlibm_TBL_atan1[index1+1]; /* lower table */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis conlo2 = __vlibm_TBL_atan1[index2+1]; /* lower table */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis poly1 = f1*((p3*tmp1 + p2)*tmp1 + p1)*tmp1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis poly2 = f2*((p3*tmp2 + p2)*tmp2 + p1)*tmp2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ansu = conup + f ; /* compute atan(f) upper */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ansu1 = conup1 + f1; /* compute atan(f) upper */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ansu2 = conup2 + f2; /* compute atan(f) upper */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ansl = (((conup - ansu) + f) + poly) + conlo ;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ansl1 = (((conup1 - ansu1) + f1) + poly1) + conlo1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ansl2 = (((conup2 - ansu2) + f2) + poly2) + conlo2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* now check to see if these are 'real' or 'dummy' arguments BEFORE storing */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *yaddr = sign ? -ans: ans; /* this one is always good */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (argcount < 3) break; /* end loop and finish up */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } while (--n > 0);