1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
/* @(#)e_acos.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25,
for performance improvement on pipelined processors.
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $";
#endif
/* __ieee754_acos(x)
* Method :
* acos(x) = pi/2 - asin(x)
* acos(-x) = pi/2 + asin(x)
* For |x|<=0.5
* acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c)
* For x>0.5
* acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
* = 2asin(sqrt((1-x)/2))
* = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z)
* = 2f + (2c + 2s*z*R(z))
* where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
* for f so that f+c ~ sqrt(z).
* For x<-0.5
* acos(x) = pi - 2asin(sqrt((1-|x|)/2))
* = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
*
* Special cases:
* if x is NaN, return x itself;
* if |x|>1, return NaN with invalid signal.
*
* Function needed: __ieee754_sqrt
*/
#include "math.h"
#include "math_private.h"
#define one qS[0]
#ifdef __STDC__
static const double
#else
static double
#endif
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
pS[] = {1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
-3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
-4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
3.47933107596021167570e-05}, /* 0x3F023DE1, 0x0DFDF709 */
qS[] ={1.0, -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
-6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
7.70381505559019352791e-02}; /* 0x3FB3B8C5, 0xB12E9282 */
#ifdef __STDC__
double __ieee754_acos(double x)
#else
double __ieee754_acos(x)
double x;
#endif
{
double z,p,q,r,w,s,c,df,p1,p2,p3,q1,q2,q3,z2,z4,z6;
int32_t hx,ix;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
if(ix>=0x3ff00000) { /* |x| >= 1 */
u_int32_t lx;
GET_LOW_WORD(lx,x);
if(((ix-0x3ff00000)|lx)==0) { /* |x|==1 */
if(hx>0) return 0.0; /* acos(1) = 0 */
else return pi+2.0*pio2_lo; /* acos(-1)= pi */
}
return (x-x)/(x-x); /* acos(|x|>1) is NaN */
}
if(ix<0x3fe00000) { /* |x| < 0.5 */
if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
z = x*x;
#ifdef DO_NOT_USE_THIS
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
#else
p1 = z*pS[0]; z2=z*z;
p2 = pS[1]+z*pS[2]; z4=z2*z2;
p3 = pS[3]+z*pS[4]; z6=z4*z2;
q1 = one+z*qS[1];
q2 = qS[2]+z*qS[3];
p = p1 + z2*p2 + z4*p3 + z6*pS[5];
q = q1 + z2*q2 + z4*qS[4];
#endif
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx<0) { /* x < -0.5 */
z = (one+x)*0.5;
#ifdef DO_NOT_USE_THIS
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
#else
p1 = z*pS[0]; z2=z*z;
p2 = pS[1]+z*pS[2]; z4=z2*z2;
p3 = pS[3]+z*pS[4]; z6=z4*z2;
q1 = one+z*qS[1];
q2 = qS[2]+z*qS[3];
p = p1 + z2*p2 + z4*p3 + z6*pS[5];
q = q1 + z2*q2 + z4*qS[4];
#endif
s = __ieee754_sqrt(z);
r = p/q;
w = r*s-pio2_lo;
return pi - 2.0*(s+w);
} else { /* x > 0.5 */
z = (one-x)*0.5;
s = __ieee754_sqrt(z);
df = s;
SET_LOW_WORD(df,0);
c = (z-df*df)/(s+df);
#ifdef DO_NOT_USE_THIS
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
#else
p1 = z*pS[0]; z2=z*z;
p2 = pS[1]+z*pS[2]; z4=z2*z2;
p3 = pS[3]+z*pS[4]; z6=z4*z2;
q1 = one+z*qS[1];
q2 = qS[2]+z*qS[3];
p = p1 + z2*p2 + z4*p3 + z6*pS[5];
q = q1 + z2*q2 + z4*qS[4];
#endif
r = p/q;
w = r*s+c;
return 2.0*(df+w);
}
}
|
