2N/A * 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. * See the License for the specific language governing permissions * and limitations under the License. * When distributing Covered Code, include this CDDL HEADER in each * 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] * Copyright 2011 Nexenta Systems, Inc. All rights reserved. * Copyright 2006 Sun Microsystems, Inc. All rights reserved. * Use is subject to license terms. #
define HI(x) *(
1+(
int*)x)
#
define LO(x) *(
unsigned*)x
#
define LO(x) *(
1+(
unsigned*)x)
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 */ 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 */ 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 */ ans = f - f;
/* return NaN if x=NaN*/ else if (
intf <
0x3e300000)
/* avoid underflow for small arg */ else if (
intf >
0x43600000)
/* avoid underflow for big arg */ *y = (
sign) ? -
ans:
ans;
/* store answer, with sign bit */ 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 */ 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) */ 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 */ f1 =
0.0;
/* put dummy values in args 1,2 */ goto UNROLL3;
/* finish up with 1 good arg */ /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ 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 */ ans =
f1 -
f1;
/* return NaN if x=NaN*/ else if (
intf <
0x3e300000)
/* avoid underflow for small arg */ else if (
intf >
0x43600000)
/* avoid underflow for big arg */ *y = (
sign1) ? -
ans:
ans;
/* store answer, with sign bit */ argcount =
1;
/* we still have 1 good arg */ f1 =
0.0;
/* put dummy values in args 1,2 */ 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 */ 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) */ 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 */ f2 =
0.0;
/* put dummy value in arg 2 */ goto UNROLL3;
/* finish up with 2 good args */ /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ 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 */ ans =
f2 -
f2;
/* return NaN if x=NaN*/ else if (
intf <
0x3e300000)
/* avoid underflow for small arg */ else if (
intf >
0x43600000)
/* avoid underflow for big arg */ *y = (
sign2) ? -
ans:
ans;
/* store answer, with sign bit */ argcount =
2;
/* we still have 2 good args */ f2 =
0.0;
/* put dummy value in arg 2 */ 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 */ 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) */ 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 ansu =
conup + f ;
/* compute atan(f) upper */ /* now check to see if these are 'real' or 'dummy' arguments BEFORE storing */ if (
argcount <
3)
break;
/* end loop and finish up */