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